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Abstract 

Explosive  Ordnance  Disposal  (EOD)  Open  Burning  (OB)  operations  are 
performed  to  treat  and  dispose  of  unserviceable  munitions  m  the  Department  of  Defense 
(DOD)  inventory.  Air  pollution  modeling  of  OB  operations  is  an  Environmental 
Protection  Agency  requirement  for  permit  issuance  at  OB  sites.  Specific  OB  regulation  is 
still  in  its  infancy;  therefore,  estabhshment  of  OB  modeling  techniques  is  still  in  the  early 
stages.  This  thesis  effort  sought  to  develop  a  computer  model,  based  upon  the  Gaussian 
Puff  Equation.  The  model  varies  firom  standard  plume  modeling  practices  by  not  making 
the  assumption  that  the  wind  direction,  wind  speed  and  turbulence  are  uniform  throughout 
the  duration  of  the  bum.  The  model  assigns  meteorological  data  to  each  explosion  (pufi) 
generated  by  the  OB  source.  The  experiments  in  this  research  effort  assigned 
meteorological  data  to  the  pirflfs  based  upon  averaging  the  weather  data  over  1,  10,  and 
60  minute  periods.  The  results  of  the  research  showed  that  there  was  a  statistically 
significant  difference  (95%  confidence)  between  1  minute  and  60  minute  weather  data 
plume  concentrations  in  the  receptor  grid  in  100%  of  the  experiments  performed.  In  each 
case,  the  one  minute  weather  data  produced  smaller  average  plume  concentrations  in  the 
receptor  grid  than  the  sixty  minute  weather  data.  This  suggests  that  a  closer  look  is 
needed  when  estabhshing  standard  OB  modeling  techniques  involving  meteorological 
input  data  to  insure  permits  are  not  denied  due  to  conservative  approaches  involving 


uniform  weather  data. 


OB  operation,  con:q)uter  modeling  is  an  acceptable  way  to  determine  plume  concentrations 
in  a  receptor  grid. 


Current  Gaussian  Puff  based  models  assume  a  steady  state  condition,  i.e.,  uniform 
wind  speed,  direction,  and  turbulence  throughout  the  time  period  of  interest.  Some  of 
these  models  also  use  a  single  stability  classification  scheme  to  determine  dispersion 
parameters  in  both  the  horizontal  and  vertical  directions. 

The  purpose  of  this  research  is  to  develop  a  FORTRAN  77  model  to  simulate  puff 
releases  from  an  EOD  open  burning  operation.  The  model  used  sets  of  meteorological 
data  averaged  over  1,10,  and  60  minutes.  Each  puff  released  fi’om  the  OB  source  was 
assigned  specific  weather  data  based  upon  the  averagiag  tune  of  interest.  The  model’s 
dispersion  parameters  were  calculated  using  separate  schemes  for  the  vertical  and 
borizontal  dispersion  parameters.  The  data  produced  by  the  model  were  used  to 
determine  whether  a  statistically  significant  difference  exists  between  the  use  of  1,  10,  and 
60  minute  meteorological  data. 

The  following  chapters  will  develop  the  argument  for  the  need  to  examine  the  use 
of  more  appropriate  averaging  times  when  modeling  EOD  Open  Burning  sources. 
Additionally,  the  need  to  examine  the  use  of  site-specific  meteorological  data  to  model 
EOD  Open  Bmning  sources  will  also  be  investigated.  Chapter  Two  reviews  the 
importance  of  OB  activities  to  DOD  and  the  necessity  to  produce  models  which  more 
closely  mimic  OB  operations.  Chapter  Three  explains  the  methodology  employed  in  this 
study  to  examine  the  effects  of  meteorological  data  averaging  times.  Chapter  Four 
presents  the  results  of  this  study.  Finally,  Chapter  Five  interprets  the  meaning  of  the 
results  and  suggests  future  research  possibilities. 
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II.  Literature  Review 


Practical  Justification 

The  stockpile  of  excess  mimitions  in  the  Department  of  Defense  (DOD)  inventory 
requiring  disposal  as  hazardous  waste  has  drastically  increased  over  the  last  ten  years.  “At 
the  end  of  1993,  DOD  estimated  the  size  of  the  munitions  stockpile  requiring  safe 
demilitarization  by  treatment  and/or  disposal  at  over  435,000  tons,  almost  doubling  its  size 
since  1982  when  the  stockpile  was  an  estimated  250,000  tons”  (Howell  and  Tope,  1994). 
“Because  of  the  enormous  backlog  of  military  weapons  that  require  timely  demditarization 
and/or  disposal .  .  .  DOD  recognizes  that  the  surplus  of  ordnance  is  fast  approaching  a 
critical  level.  .  .”  (Tope  and  Howell,  1994).  There  are  two  main  reasons  for  this  marked 
increase  in  the  size  of  this  stockpile  of  unserviceable  mimitions.  First  is  the  recent 
downsizing  of  mditary  operations  following  the  end  of  the  Cold  War.  Downsizing  has 
reduced  the  number  of  military  installations,  operations,  and  weapons  systems  which  hi 
turn  increased  the  rate  of  unserviceable  munitions  generated  (Tope  and  Howell,  1994). 
Second,  there  was  and  still  is  a  great  deal  of  confusion  between  the  regulatory  agencies 
and  the  regulated  community  regarding  classification  schemes  for  munitions.  Siace  it 
appears  the  mihtary  is  still  in  a  downsizmg  mode,  that  source  of  increase  of  the  munitions 
stockpile  is  not  going  away  in  the  near  future.  Therefore,  the  regulatory  confusion  will  be 
investigated  as  a  way  to  decrease  the  stockpile.  It  is  necessary  to  address  exactly  how 
munitions  are  handled  and  ultimately  destroyed  before  the  air  quahty  issues  associated 
with  their  destruction  are  discussed. 

The  Resource  Conservation  and  Recovery  Act  (RCRA)  governs  the  regulation  of 
hazardous  waste.  Once  an  item  can  no  longer  be  used  for  its  intended  purpose,  it  is 
subject  to  RCRA  which  identifies  two  classification  schemes  for  determining  whether  an 
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item  is  hazardous  waste.  Either  the  waste  is  hsted,  i.e.,  the  item  is  specifically  mentioned 
on  an  Environmental  Protection  Agency(EPA)  list,  or  it  is  a  characteristic  waste. 
Characteristic  wastes  display  any  the  following  properties:  toxicity,  ignitahihty, 
flammabihty  or  reactivity.  Once  classified  as  hazardous  waste,  the  item  is  subject  to 
stringent  requirements,  including  a  storage  time  limit  of  90  days  without  a  special  permit. 
Herein  hes  the  crux  of  the  problem  with  the  munitions  stockpile  --  storage  time. 

Mditary  munitions  “generally  consist  of  an  assortment  of  explosive  [i.e.,  reactive] 
fin  materials,  associated  metal  and  plastic  casings,  projectiles  and  primer  components” 
(Tope  and  Howell,  1994).  Common  examples  of  munitions  include:  bullets,  bombs, 
ejection  seat  cartridges,  flares  and  smoke  grenades.  The  base  Munitions  Accountable 
Systems  Officer  (MASO)  closely  monitors  the  munitions  inventory  by  accomphshing  daily 
inventory  inspections  to  determine  the  condition  of  each  item  Condition  codes  are 
assigned  to  each  item  based  upon  whether  the  item  is  obsolete,  excess  (expired  shelf-life), 
or  damaged  (leaking,  missing  parts,  etc.)  (Johnson,  1994).  Once  a  MASO  declares  an 
item  to  be  unserviceable  for  any  reason,  the  process  of  deciding  how  to  ultimately  deal 
with  the  item  begins. 

There  are  three  primary  options  for  dealing  with  unserviceable  munitions.  First, 
munitions  supply  personnel  determine  whether  the  item  can  be  used  at  another  installation. 
The  second  option  is  the  use  of  the  item  for  legitimate  training  by  mihtary  personnel.  The 
third  option  (and  last  resort)  is  a  condition  code  H  designation.  Condition  code  H 
signifies  the  item  is  ‘unserviceable’  and  therefore  must  be  designated  as  RCRA  hazardous 
waste  since  it  can  no  longer  be  used  for  its  intended  purpose.  Condition  code  H  items  are 
physically  separated  from  serviceable  munitions  and  the  MASO  schedules  the  items  for 
disposal  (HQ  USAF  Message,  1994).  The  ‘RCRA  90-day  storage  clock  starts  ticking’  as 
soon  as  the  item  is  declared  condition  code  H.  Within  the  90-day  time  limit,  personnel 
must  complete  volumes  of  mihtary  paperwork  intended  to  control  these  sensitive. 
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dangerous  items,  arrange  for  transportation  to  an  authorized  disposal  site,  and  actually 
transport  the  items  to  the  disposal  site.  Ultimately,  RCRA  requires  mrmitions  to  be 
disposed  of  by  means  of  open  burning  (OB)  or  open  detonation  (OD)  within  90  days  of 
the  initial  designation  as  condition  code  H,  i.e.  hazardous  waste. 

In  accordance  with  Air  Force  regulations,  highly  trained  Explosive  Ordnance 
Disposal  (EOD)  personnel  determine  whether  the  item  will  be  destroyed  by  OB  or  OD. 
EOD  personnel  must  follow  strict  regulations  which  specifically  outline  the  proper  method 
of  disposal,  i.e.  open  burning  or  open  detonation,  for  each  item  in  the  Air  Force  munitions 
inventory.  Due  to  the  fact  that  OB/OD  operations  involve  the  disposal  of  hazardous 
waste,  the  installations  where  OB/OD  is  conducted  are  regulated  under  RCRA  as 
treatment,  storage,  and  disposal  facilities  (TSDFs).  This  study  focuses  on  destruction  of 
munitions  by  open  bunung. 

Very  simphstic  ia  nature,  the  OB  process  evolved  over  the  past  40  years  into  the 
preferred  munitions  destruction  method  for  items  that  cannot  be  destroyed  by  open 
detonation.  EOD  personnel  begin  OB  events  by  stacking  wooden  cargo  pallets  in  a 
shallow  pit,  bermed  area  or  bum  pan  in  order  to  provide  fuel  to  keep  the  fire  burning. 

Next,  EOD  piles  the  items  to  be  destroyed  on  top  of  the  dunnage  and  pours  either  jet  or 
diesel  fuel  on  top  of  the  pile.  Finally,  they  hght  the  pile  using  initiator  charges  timed  to 
allow  EOD  personnel  to  clear  the  area. 

The  explosions  involved  in  OB  activities  are  different  fi'om  conventional  ideas 
concemmg  explosions  in  one  significant  manner.  An  OB  event  does  not  involve  one  large 
explosion.  Instead,  OB  events  involve  a  series  of  small  explosions  as  the  munitions  items 
being  destroyed  reach  temperatures  high  enough  to  produce  the  exothermic  (heat  and 
energy  releasing)  reaction  needed  to  thoroughly  destroy  the  items  (Tope  and  HoweU, 
1994). 


2-3 


OB  operations  are  only  one  aspect  of  the  mission  of  EOD  personnel.  Other  duties 
include  providing  security  for  VIPs,  emergency  response  to  bomb  threats,  training,  and 
classified  missions.  These  other  duties  occupy  a  great  deal  of  the  EOD  personnel’s  time; 
therefore,  OB  operations  only  occur  on  a  monthly  or  quarterly  basis.  The  rarity  of  OB 
events  increases  the  size  of  the  munitions  stockpile  because  the  supply  cannot  keep  up 
with  the  demand,  i.e.,  there  are  not  enough  OB  facilities  and  EOD  personnel  to  meet  the 
disposal  demand.  Therefore,  in  order  to  reduce  the  stockpile  and  keep  it  at  a  manageable 
level,  it  is  imperative  that  more  facUities  be  fully  permitted  in  a  timely  manner. 

EPA  RCRA  permitting  is  a  long  process  due  to  the  large  number  of  apphcations 
received.  The  Department  of  Defense’s  experience  with  the  EPA’s  OB  permit  apphcation 
process  began  in  1980  and  to  date,  only  one  final  permit  has  been  issued  (TNRCC  and 
EPA,  Apr  and  June  1992).  The  delay  in  permitting  of  OB  units  is  due  in  part  to  the 
uniqueness  of  this  disposal  method. 

OB  units  are  regulated  in  the  CFR  as  ‘Miscellaneous  Subpart  X  Units’  (EPA, 
1987).  The  CFR  specifies  design  and  operating  standards  for  more  common  hazardous 
waste  disposal  methods,  while  standards  for  Subpart  X  units,  such  as  OB  units,  are  permit 
specific.  The  EPA’s  lack  of  sufficient  knowledge  to  promulgate  specific  design  and 
operating  standards  for  OB  causes  major  delays  in  the  permitting  process  (Federal 
Register,  1987;  Johnson,  1994).  Except  for  the  single  permitted  facility,  other  OB  units 
operate  under  temporary  ‘interim-status’  permits  until  final  permits  are  issued  or  the 
interim- status  expires  and  the  unit  is  set  up  for  closure  (Davis,  1995). 

Initially,  a  great  deal  of  confusion  existed  regarding  several  issues  surrmmding  OB 
operations.  This  confusion  prompted  several  concerned  wing  commanders  to  temporarily 
shut  down  their  OB/OD  units  until  definitive  guidance  was  formulated  by  DOD. 

Justifiably  concerned,  these  commanders  thought  this  regulatory  confusion  might  lead  to 
the  receipt  of  a  Notice  of  Violation,  which  could  be  accompanied  by  a  fine,  from 
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eaviromnental  regulators.  Tlie  temporary  shutdown  of  these  units  caused  a  bottleneck  in 
the  unserviceable  munitions  disposal  process  resulting  in  even  more  condition  code  H 
munitions  items  awaiting  disposal. 

Currently,  RCRA  and  the  Clean  Air  Act  do  not  require  a  separate  air  permit  to 
operate  an  OB  thermal  treatment  unit.  RCRA  permits  incorporate  air  emissions 
regulations  by  using  umbrella  statements  such  as,  ‘permittee  shall  con:q)ly  with  all 
applicable  federal,  state  and  local  environmental  regulations.’  This  all  encompassing 
statement  does  not  specifically  mention  air  quahty  aspects  such  as  limits  on  certain 
chemical  constituents  in  the  smoke  plume  generated  by  OB.  Therefore,  it  would  be 
difficult  to  pick  out  this  requirements  while  wading  through  the  multitude  of  enumerated 
hazardous  waste  requirements  for  comphance  with  the  permit. 

Due  to  the  obscurity  of  air  quality  issues  in  the  past,  the  issue  of  air  emissions 
concerns  has  not  been  addressed  until  recently  when  EOD  operations  advanced  to  the 
forefront  in  the  regulatory  community.  The  growing  stockpile  of  unserviceable  munitions 
and  the  fifteen-plus  years  already  spent  toward  applying  for  OB  permits  necessitates 
attention  to  and  compliance  with  aU  environmental  aspects  of  OB,  including  air  quality 
standards. 

Unfortunately,  since  the  central  focus  of  the  regulations  governing  OB  operations 
and  mxmitions  inventory  control  centers  around  RCRA,  the  issue  of  the  air  emissions 
associated  with  burning  operations  is  sometimes  hidden  from  view.  The  new  attention 
paid  to  OB  operations  has  made  the  issue  of  air  pollution  modeling  a  permitting 
requirement.  Therefore,  it  is  necessary  to  assess  the  compliance  of  OB  units  with  air 
quality  standards  through  the  use  of  research  methods. 
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Academic  Justification 

Air  sampling  during  an  OB  operation  is  difi&cult  at  best.  The  explosive  nature  of 
the  material  being  burned  necessitates  that  personnel  in  the  area  stay  a  safe  distance  away. 
During  the  bum  operation,  certain  items  may  ‘kick  out’  of  the  burning  pile  creating  a 
dangerous  situation  for  anyone  caught  too  close  to  the  bum  area.  Unfortunately,  these 
kick-outs  also  make  the  use  of  air  sampling  equipment  prohibitive,  since  sample  equipment 
could  be  destroyed  by  a  projectile.  Consequently,  alternate  methodologies  must  be 
established  for  determining  comphance  with  air  quality  standards.  This  study  focuses  on 
modeling  plume  concentrations  generated  in  a  receptor  grid  from  an  OB  source. 

Since  first  being  suggested  in  detail  by  Pasquih  in  1961,  air  pollution  modeling 
using  the  Gaussian  distribution  has  developed  well-estabhshed  roots  (Bowen,  1994).  The 
Gaussian  Puff  Model  (Eqn  1)  is  used  to  determine  plume  concentrations  at  a  given 
receptor  location  coordinate  at  a  given  time  (x,  y,  z,  t).  This  model  is  used  for  cases 
involving  instantaneous  releases,  i.e.,  an  explosion  (Turner,  1969). 

Currently,  there  is  no  EPA  recommended  model  to  handle  the  special  set  of 
circumstances  surrounding  OB  sources.  In  1994,  a  model  development  program  was 
undertaken  by  the  Department  of  Defense  (DOD)  and  the  Department  of  Energy  (DOE) 
to  develop  a  model  capable  of  handling  OB  sources  (Weil  and  Templeman,  1995). 

Standard  models  are  based  upon  a  Gaussian  Puff  Model  (Eqn  1)  which  assumes 
uniform  average  ■wind  speed  (  m  ),  wind  direction  and  turbulence  throughout  the  duration 
of  the  OB  operation  (Bowen,  1994;  GrifBn,  1994).  The  meteorological  input  data  for 
these  models  is  usually  measured  at  a  height  of  10m  over  an  averaging  time  of  one  hour 
(but  not  less  than  10  minutes)  (MitcheU,  1982).  These  models  attempt  to  account  for  the 
changes  experienced  by  the  atmosphere  through  the  use  of  the  dispersion  parameters  to 
describe  material  spread  in  the  puffs.  This  crude  method  does  not  sufficiently  account  for 
the  constantly  changing  wind  speed,  direction  and  turbulence  experienced  by  our 
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atmosphere.  Examination  of  Equation  1  shows  that  these  atmospheric  characteristics  have 
a  direct  effect  on  modeled  plume  concentrations  in  a  receptor  grid. 


C{x,y,z,t)  = 


207’ 


{2n) 


3/2 


-exp 


^  z 


^ -{x  -  u 


2(7 


y 


2ct 


2c7lj 


(1) 


The  source  term  (Qt)  in  the  Gaussian  Puff  Model  is  the  total  mass  of  the  material 
released  during  the  time  period  of  interest  (Turner,  1969).  This  type  of  source  term  is 
vahd  when  one  explosion  is  involved,  but  would  appear  to  he  invalid  when  applying  the 
puff  model  to  OB  operations  involving  a  series  of  explosions. 

The  horizontal  (ctx  and  oy)  and  vertical  (oz)  dispersion  parameters  in  the  standard 
Gaussian  Puff  Model  are  uniform  throughout  the  time  period  of  interest  (usually  1  hour). 
Furthermore,  typical  puff  models  use  a  single  dispersion  parameter  classification  scheme  to 
determine  both  the  horizontal  and  vertical  parameters.  The  use  of  uniform  dispersion 
parameters  throughout  the  time  period  of  interest  may  he  appropriate  in  situations 
involving  one  explosion,  but  is  clearly  inappropriate  for  OB  activities  involving  a  series  of 
explosions. 

Several  schemes  exist  for  the  determination  of  the  dispersion  parameters  (ox,  Oy, 
and  Gz).  The  most  frequently  used  schemes  are  the  Vertical  Temperature  Gradient,  Sigma 
Theta,  and  Modified  Sigma  Theta.  The  Vertical  Temperature  Gradient  (DT/DZ)  involves 
the  use  of  the  difference  in  temperatures  between  two  specified  to  determine  the  stabdity 
conditions  and,  subsequently,  the  dispersion  parameters.  The  Sigma  Theta  method  uses 
the  standard  deviation  of  the  horizontal  wind  direction  to  determine  stabdity  conditions. 
The  Modified  Sigma  Theta  (MST)  also  uses  the  standard  deviation  of  the  wind  direction 
to  determine  stabdity  conditions.  However,  MST  accounts  for  nighttime  conditions  by 
factoring  hi  the  wind  speed  when  determining  the  stabdity  class  from  which  to  calculate 
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the  dispersion  parameters.  The  MST  is  inappropriate  for  use  in  modeling  OB  activities 
since  safety  considerations  dictate  that  OB  activities  occur  only  during  dayhght  hours. 

Models  developed  for  air  pollution  typically  use  one  atmospheric  stabihty 
classification  method  to  calculate  all  three  stabihty  parameters.  This  means  that  if  DT/DZ 
is  used,  assumptions  about  horizontal  stabihty  class  would  be  made  based  upon 
measurements  taken  in  the  vertical  direction,  i.e.,  temperature  gradient.  Additionahy,  if 
Sigma  Theta  is  used,  assumptions  about  vertical  stabihty  class  would  be  made  based  upon 
measurements  taken  in  the  horizontal  direction,  i.e.  wind  direction.  Puff  models  which  do 
not  assume  (cJx  =  Oy  -  Cz)  wih  often  assume  that  the  puffs  released  from  the  source  have 
a  circular  horizontal  cross  section;  therefore,  (CJx  =  CJy). 

Modification  of  horizontal  dispersion  parameters  due  to  an  averaging  time  different 
than  the  ten  minutes  estabhshed  by  the  Pasquih-Gifford  curves  are  typicaUy  modified  by 
the  one-fifth  power  law.  The  one- fifth  power  law  is  used  to  modify  data  from  ten  minute 
averaging  time  to  another  averaging  time  of  interest  (Kunkel,  1991).  This  law  was 
derived  from  empirical  data  and  may  not  be  vahd  in  ah  situations.  Obtahung  the 
meteorological  data  needed  to  calculate  the  dispersion  parameters  is  preferable  to  using  an 
estimated  formula  to  calculate  the  dispersion  parameters  for  differing  averaging  times. 

This  formula  is  necessary  when  the  meteorological  data  set  is  not  available  to  perform 
calculations  for  differing  averaging  times. 

The  availabihty  of  site-specific  meteorological  weather  data  makes  the  modeling  of 
OB  operations  more  reahstic.  The  creation  of  mobile  meteorological  towers  enables 
researchers  to  gather  site-specific  data  versus  the  common  practice  of  calling  the  local 
airport  to  obtain  weather  data  (Weil,  1995).  The  weather  data  received  from  the  local 
airport  may  be  vahd  for  the  immediate  area  around  the  airport,  but  may  be  vastly  different 
from  the  weather  at  the  location  of  the  source. 
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EPA  requirements  for  EOD  Open  Burning  permit  applications  are  likely  to  include 
site-specific  assessments  which  assign  risk  numbers  to  facilities.  The  varying  locations  of 
OB  sources  necessitates  the  use  of  site  meteorological  data  to  insure  the  most  accurate 
modeling  is  accomphshed  for  that  site. 
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lil.  Methodology 


FORTRAN  Computer  Model 

Due  to  the  inherent  difficulty  and  danger  involved  in  performing  air  measurements 
above  an  Explosive  Ordnance  Disposal  (EOD)  Open  Burning  (OB)  operation, 
mathematical  modeling  must  be  used  to  predict  plume  concentrations.  A  modified 
Gaussian  Puff  Equation  was  used  to  develop  a  FORTRAN  computer  model  named 
PUFFY  (See  Appendix  A  for  FORTRAN  code).  PUFFY  is  used  to  calculate  pollutant 
concentrations  in  a  Cartesian  receptor  grid  extending  1  km  from  the  source  (the  OB 
operation)  hi  all  four  Cardinal  directions,  North,  South,  East,  and  West.  PUFFY  was 
designed  with  the  assertion  that  the  individual  pufis  released  during  the  eiqplosions 
occurring  during  an  OB  operation  may  potentially  travel  along  various  trajectories  (Figure 
1)  versus  the  traditional  assunq)tion  that  puffs  would  travel  along  the  same  trajectory. 


Z 


3-1 


Origin  of  Meteorological  Data 

The  meteorological  data  used  as  inputs  for  the  model  developed  in  this  thesis  eflfort 
was  gathered  from  five  weather  towers  at  various  locations  in  and  arormd  Cedar  Bog  in 
(Springfield,  OH)  as  part  of  an  on-going  research  effort  by  Wright  State  University 
(WSU).  The  towers,  purchased  from  Clrmatronics,  Inc.,  measiue  the  following 
parameters  relevant  to  this  thesis  effort;  temperature,  wind  speed,  wind  direction  and 
standard  deviation  of  the  wind  direction  (00)  at  Im  and  10m  heights.  The  towers  were  set 
using  a  thermosistor  to  measure  temperature,  anemometers  to  measure  wind  speed  and  a 
horizontal  wind  vane  to  measure  the  wind  directions  used  to  calculate  the  standard 
deviation.  Originally  set  to  take  weather  data  at  15  minute  intervals,  the  tower 
measurement  interval  was  reset  to  one  minute  to  gather  data  for  this  thesis  effort.  The  one 
minute  data  was  taken  for  a  period  of  two  days. 

The  data  gathered  from  the  Climatronics  towers  was  manipulated  by  a  FORTRAN 
program  called  ‘METALL’  so  that  it  would  be  in  a  usable  form  for  the  FORTRAN 
program  developed  for  this  thesis  effort.  The  program  METALL  was  used  to  convert  the 
relevant  weather  data  from  1  minute  averaging  time  to  10  minute  and  60  minute  averages. 

The  program  separates  the  wind  directions  into  x  and  y  conqjonents  during  the 
time  period  of  interest,  i.e.,  ten  minutes  or  sixty  minutes,  and  totals  them  Separate  x  and 
y  component  averages  are  calculated  and  the  resulting  wind  direction  is  computed  by 
taking  the  inverse  tangent  of  the  result  of  the  y  component  divided  by  the  x  component. 
Next,  180  degrees  was  then  added  to  the  resulting  angle  due  to  meteorological  convention 
for  wind  direction.  Meteorological  convention  describes  the  direction  from  which  the 
wind  originated,  i.e.,  from  the  South.  Wind  from  due  South  is  considered  the  0  degree 
starting  point  with  wind  directions  measured  in  a  clockwise  direction.  In  contrast, 

Cartesian  coordinates  describe  wind  direction  in  terms  of  the  direction  to  which  it  is 
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blowing,  i.e.,  wind  blowing  toward  the  East  is  considered  the  0  degree  starting  point. 
Successive  wind  directions  are  measured  in  a  counterclockwise  direction. 

Instead  of  using  the  weather  tower  values  of  09  directly  to  eventually  determine 
the  value  of  (ay)c  by  using  the  one  fifth  power-law  (Eqn  1)  to  convert  the  Oy  values  to  the 
appropriate  averaging  times,  the  derivation  in  Appendix  C  was  used  to  justify  the  formula 
(Eqn  2)  used  to  calculate  the  (09)  for  the  appropriate  averaging  times  based  upon  one 
minute  averaging  time  weather  tower  data. 


(a  y)^=  a  y{tc  n^f'^  (1) 

where, 

(cfY)c  =  Corrected  Horizontal  Dispersion  Parameter  (meters) 
tc  =  Weather  Data  Averaging  Time  (minutes) 


(2) 


where, 

Nx  =  Ratio  of  the  number  of  readings  taken  in  one  minute  to  the  number  of  readings 
taken  in  the  desired  averagmg  time  (i.e.,  for  10  minutes— 60/600=1/10) 

09c  =  Corrected  Standard  Deviation 

09  =  Standard  Deviation  of  the  Wind 

p.  =  Individual  Wind  Directions  (from  Climatronics  Towers) 

yU  =  Average  Wind  Direction  for  Time  Period  of  Interest  (10  or  60  minutes) 


Source  Term 

Each  experiment  conducted  used  a  one  hour  time  period  of  interest.  The  Open 
Bum  activity  was  simulated  by  the  release  of  one  puff  during  each  minute.  Each  puff  is  of 
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uniform  strength.  The  model  has  the  built-in  flexibihty  for  the  user  to  change  the  strength 
of  the  puffs  released.  The  user  also  has  the  capabihty  to  change  the  source  strength  so 
that  each  puff  released  has  a  different  strength. 

Calculation  of  Diffusion  Parameters 

Horizontal  Diffusion  Parameter.  The  standard  deviation  of  the  horizontal  wind 
direction  (oq)  was  used  to  determine  the  stabihty  parameter  needed  to  calculate  the 
horizontal  diffiision  parameters  and  Oy).  The  Turner  Workbook  describes  more 
detailed  information  about  the  concept  of  the  diffiision  parameter  (Turner,  1969).  First, 
the  value  of  (sigma  theta)  is  read  in  from  the  data  file  generated  by  the  FORTRAN 
program  ‘METAfLL.’  Sigma  Theta  (a©)  is  then  used  to  determine  the  value  of  the  stabihty 
parameter.  Table  1  illustrates  the  discrete  values  commonly  used  to  determine  stabihty 
parameters  fi'om  ranges  of  ae  values. 


TABLE  1  PasquUl  Stabihty  Categories  and  Corresponding  Stability  Parameter  Values  with 
Associated  Ranges  of  (Sigma  Theta)  --  Using  the  Modified  Sigma  Theta  (MST)  Method 


A 

B 

C 

D 

E 

F 

Stabihty  Parameter 

0.5 

1.5 

2.5 

3.5 

4.5 

5.5 

Sigma  Theta 

>22.5 

17.5  to 

12.5  to 

7.5  to 

3.8  to 

<3.8 

(degrees) 

22.5 

17.5 

12.5 

7.5 

This  thesis  effort  asserts  that  since  the  atmosphere  does  not  normaUy  behave  in  a  discrete 
manner,  the  values  of  the  stability  parameters  used  in  the  puff  model  should  not  be 
discrete.  Linear  interpolation  between  stabihty  parameters(SP)  was  used  to  determine  a 
more  precise  value  of  the  stabihty  parameter  than  the  estabhshed  practice  of  choosing 
discrete  stabihty  classes  and  a  corresponding  discrete  stabihty  parameters.  The  value  of 
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ay  is  determined  by  the  following  power  law  expression  (Equation  3)  derived  from  the 
Pasquill- Gifford  curves  (Figure  2)  (Seinfeld,  1986;  Kunkel  1991) 

cr^  =  ax^  (3) 

where, 

a  =  0.479  -  0. 1232*SP  +  0.00904*SP^ 

SP  =  Stability  Parameter 

b  =0.9 

The  assumption  was  made  that  the  puff  released  from  the  EOD  operations  has  a  circular 
horizontal  cross-section;  therefore,  <y  ^  y 
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FIGURE  2  Correlations  for  Cy  based  on  the  Pasquill  stabUity  classes  A-F  (Seinfeld,  1986) 
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Vertical  Diffusion  Parameter.  The  vertical  temperature  gradient  (dT/dZ)  was 
used  to  determine  the  stability  parameter  needed  to  calculate  the  vertical  dijSusion 
parameter  (az).  First,  the  temperatures  at  the  10m  and  Im  heights  were  converted  from 
Fahrenheit  to  Celsius.  Second,  the  temperature  at  Im  was  subtracted  from  the 
temperature  at  10m  (dT).  Next,  dT  was  divided  by  the  difference  in  heights  where  the 
temperatures  were  measured  (dZ)  to  determine  dT/dZ. 

Table  2  illustrates  the  values  for  the  stabihty  parameter  assigned  by  F.  V.  Hansen 
to  each  of  the  Pasquill  stabihty  categories  (Kimkel,  1991).  Table  2  also  depicts  the  ranges 
of  temperature  gradient  values  associated  with  each  stabihty  category. 


TABLE  2  Stabihty  Parameter  Values  Associated  with  Temperature  Gradient  Ranges  and 

Corresponding  Pasquill  Stabihty  Categories 


A 

B 

C 

D 

E 

F 

Stabihty  Parameter 

0.5 

1.5 

2.5 

3.5 

■9 

5.5 

DT/DZ 

<-1.9 

-1.9  to 

-1.7  to 

-1.5  to 

-0.5  to 

>+1.5 

(degrees  C/lOOm) 

-1.7 

i._  -Lj 

-0.5 

+1.5 

Once  the  value  of  DT/DZ  was  calculated,  it  was  used  to  determine  the  value  of  the 
stabihty  parameter  based  upon  the  ranges  of  the  temperature  gradient  recommended  for 
use  by  the  Nuclear  Regulatory  Commission  (NRC)  method  for  determining  vertical 
stabihty  class.  Using  the  value  of  the  calcidated  stabihty  parameter,  the  values  of  the 
coefficient,  c,  and  exponent,  d,  are  computed.  The  values  of  the  coefficient  and  exponent 
are  needed  to  compute  az  according  to  the  power-law  expression  in  Equation  4. 
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(4) 


The  values  of  the  coefficient  and  exponent  for  the  power-law  expression  were  determined 
by  the  use  of  the  discrete  values  in  Table  3.  The  values  of  the  stabUity  parameters  were 
reconfigiued  as  ranges  in  order  to  account  for  the  continuous  nature  of  the  stabihty 
parameters  generated  by  the  use  of  the  temperature  gradient  values. 


TABLE  3  Coefficients  and  Exponents  for  the  Vertical  Dispersion  Parameter 


A 

B 

C 

D 

E 

F 

Stabihty  Parameter 

<1.0 

1.0  to 

2.0  to 

3.0  to 

4.0  to 

>5.0 

2.0 

3.0 

4.0 

5.0 

wmm 

745 

745 

2000 

1100 

1400 

1400 

c 

0.0414 

0.1036 

0.1173 

0.0975 

0.1050 

0.0617 

d 

1.3155 

1.0026 

0.9112 

0.8414 

0.7692 

0.7884 

x(m)> 

745 

745 

2000 

1100 

1400 

1400 

c 

1.928E-04 

0.0534 

0.4422 

0.6097 

0.8788 

0.9990 

d 

2.1234 

1.1029 

0.7382 

0.5808 

0.4771 

0.4771 

Unlike  the  coefficient,  a,  for  the  horizontal  diffiision  parameters,  a  reasonable  line  of  best 
fit  could  not  be  found  for  the  coefficient,  c,  used  to  calculate  the  vertical  dffiusion 
parameter  (Figure  3).  The  curves  used  to  calculate  the  coefficient,  c,  are  not  straight  lines 
even  when  plotted  on  logarithmic  paper.  Therefore,  linear  interpolation  between  the 
values  of  the  coefficient,  c,  and  the  exponent,  d,  would  be  an  inappropriate  approach. 
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FIGURE  3  Correlations  for  Gz  based  on  the  Pasquill  stability  classes  A-F  (Seinfeld,  1986) 

Assignment  of  Meteorological  Data  to  Individual  Puffs 

Each  puff  released  during  the  time  period  of  interest  was  assigned  its  own  set  of 
meteorological  data  depending  upon  the  averaging  time  under  study.  If  a  puff  is  released 
with  one  minute  weather  data,  each  puff  is  assigned  a  different  set  of  one  minute  weather 
data,  i.e.,  the  weather  changes  every  minute.  If  a  puff  is  released  with  ten  minute  weather 
data,  each  set  often  puffs  released  during  the  time  period  of  interest  will  have  a  separate 
set  often  minute  weather  data,  i.e.,  the  weather  only  changes  every  ten  minutes.  If  the 
puff  is  released  with  sbcty  minute  weather,  each  puff  will  be  assigned  the  same  weather 
data,  i.e.  the  weather  does  not  change  in  a  one-hour  period 
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X  and  Y  Coordinate  Determination 

The  Gaussian  Instantaneous  Plume  Equation  assumes  that  the  wind  is  blowing 
along  the  x-axis  of  a  ‘right-hand’  three-dimensional  coordinate  system.  The  coordinates 
of  the  receptor  point  of  interest  (x,  y,  z)  represent  the  distance  downwind  of  the  source(x), 
the  distance  from  the  centerline  of  the  puff  to  the  receptor  (y),  and  the  vertical  distance 
from  ground-level  to  the  receptor  (z). 

The  computer  used  to  calculate  the  plume  concentrations  in  the  PUFFY  model 
does  not  use  the  meteorological  convention  for  wind  directions,  i.e.,  in  terms  of  where 
they  origmated.  Instead,  it  uses  the  Cartesian  coordinate  system  by  manipulatmg  angles  in 
terms  of  which  direction  the  wind  is  blowing.  In  order  to  insure  the  proper  coordmates 
were  used  to  determine  the  x,  y,  and  z  distances  from  the  source  to  the  receptor,  the  wmd 
direction  from  the  weather  data  set  was  converted  to  Cartesian  coordinates.  The 
transformation  was  accomphshed  by  first  adding  180  degrees  to  the  wind  direction  of 
interest  in  order  to  represent  the  wmd  direction  as  blowing  toward  a  direction  instead  of 
resulting  from  a  direction.  Second,  the  direction  was  separated  into  its  x  and  y 
components.  Third,  the  x  component  and  y  component  were  reversed  to  account  for  the 
meteorological  convention  of  measuring  angles  in  a  clockwise  rather  than 
counterclockwise  (Cartesian)  manner.  Finally,  the  inverse  tangent  fimction  was  used  to 
determine  the  wind  direction  (Alpha)  in  Cartesian  coordinates. 

Since  the  PUFFY  model  assigns  each  puff  released  to  a  set  of  meteorological  data, 
it  was  necessary  to  consider  the  receptor  grid  used  in  the  experiment  as  a  stationary 
(reference)  grid  which  did  not  change  when  a  new  set  of  weather  data  was  initiated. 
Instead,  with  the  release  of  each  puff  the  coordinates  in  the  receptor  grid  were 
transformed  into  coordmates  for  the  wind  direction  assigned  to  the  puff  of  interest.  The 
angle  (Beta)  between  the  reference  frame’s  x-axis  and  the  receptor  grid  point  was 
determined  in  order  to  calculate  the  transformed  coordmates.  The  difference  (Gamma) 


3-9 


between  the  angles  (Alpha  and  Beta)  was  used  to  calculate  the  transformed  coordinates. 
The  transformed  coordinates  were  used  as  the  x  and  y  values  needed  to  calculate  the 
plume  concentrations. 

Reference  Grid  Spacing  and  Time  Step.  The  receptor  grid  used  for  this  thesis 
effort  consisted  of  grid  points  every  50  m  extending  out  to  1  km  from  the  Explosive 
Ordnance  Disposal  (EOD)  Open  Bmuing  operation  in  all  four  Cardinal  directions,  i.e., 
North,  South,  East,  and  West,  making  the  grid  (1  km  x  1  km).  The  grid  spacing  was  set  at 
50m  to  ensure  the  puffs  released  from  the  source  were  over  a  grid  point  in  close  pro?dmity 
to  the  source.  If  the  grid  spacing  was  any  larger,  i.e.,  100m,  a  puff  release  may  not  have 
been  detected  until  the  puff  was  well  downwind. 

The  time  step  was  set  at  15  second  intervals  to  insure  that  puffs  do  not  pass  over 
the  receptor  site  between  measurements.  For  example,  if  a  puff  has  not  yet  reached  the 
receptor  grid  point  of  interest,  at  the  10  minute  point  the  plume  concentration  is  zero.  If 
the  puff  passes  by  the  grid  point,  then  when  a  measurement  is  taken  at  the  11  minute 
point,  the  plume  concentration  is  still  zero.  This  coidd  lead  to  serious  error  if  these 
measmements  were  used  to  determine  receptor  dosage  or  maximum  exposure 
concentrations. 

Gaussian  Instantaneous  (Puff)  Equation* 

A  modification  of  the  Gaussian  Instantaneous  (Puf^  equation  was  used  to 
calculate  the  plume  concentrations  in  the  receptor  grid.  This  equation  is  often  used  in 
cases  involving  an  explosion.  The  ‘Puff  Equation’  assumes  the  source  term,  Qt, 
represents  the  total  mass  of  material  released  during  the  time  period  of  interest  (Turner, 
41).  Equation  5  is  the  Gaussian  Puff  Equation  to  calculate  plume  concentrations  from 
instantaneous  sources. 
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where, 

Qt  =  Total  mass  of  material  released  in  the  time  period  of  interest  (mass)* 

K  =3.141592654 

=  Diffusion  parameter-standard  deviation  of  material  concentration  in  the  X 
direction 

==Dij3usion  parameter- Standard  deviation  of  material  concentration  in  the  Y 
direction 

o'z  =  DiflEusion  parameter- Standard  deviation  of  material  concentration  in  the  Z 

direction 

X  =  Distance  downwind  from  the  source 

y  =  Crosswind  distance  from  the  centerline  of  the  plume  trajectory 
z  =  Vertical  distance  from  the  ground  to  the  receptor 

u  =  Wind  speed  of  individual  puff  at  the  time  of  release  (distance/time)*  at  10m 
t  =  Time  since  puff  was  released  from  source 

H  =  Effective  stack  height  (includes  the  actual  stack  height  +  plume  rise  from 
buoyant  forces) 

**  The  puff  is  assumed  to  have  a  circular  horizontal  cross  section;  therefore,  cr  x  =  y 


Statistical  Analysis 

The  model,  PUFFY,  developed  for  this  thesis  effort  incorporates  a  modified  version  of  the 
Gaussian  Puff  Equation.  Instead  of  assuming  aU  material  is  released  at  once,  PUFFY 
assumes  that  each  puff  released  has  a  separate  contribution  to  the  concentration 
measurements  at  each  grid  receptor  point  depending  upon  the  meteorological  conditions 
when  the  puff  was  released.  PUFFY  assigns  separate  meteorological  data  to  each  puff. 
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Depending  upon  the  averaging  time  of  interest,  for  one  minute  averaging  time,  unique 
meteorological  data  is  assigned  to  each  puff  during  the  one  hour  time  period  of  interest 
used  in  this  experiment.  For  ten  minute  averaging  time,  the  first  set  of  ten  puffs  have  the 
same  meteorological  data,  the  second  set  of  ten  puffs  have  a  different  set  of  data,  and  so 
on  until  the  sixth  set  of  puffs  has  been  released.  When  the  averaging  time  is  sixty  minutes, 
each  of  the  puffs  released  have  the  same  meteorological  data.  Equation  6  was  used  to 
determme  the  maximum  cumulative  concentration  to  occur  at  the  grid  point  (x,  y,  z=0). 


C(x,>;,0,0=  'Z  7^ 

p=l  (2^)- 


-{x-Uptpf 


^  X  n  ^  V  ^  rt 

yp  *'P 


where, 


p  =  Conner  for  order  in  which  puffs  were  released 

Qp  =  Total  mass  of  material  released  in  the  time  period  of  interest  (mass/time)* 

n  =3.141592654 

^xp  =  Diffusion  parameter- Standard  deviation  of  material  concentration  in  the  X 
direction  for  puff  of  interest 

=Diffusion  parameter- Standard  deviation  of  material  concentration  m  the  Y 
direction  for  puff  of  interest 

^zp  =  Diffusion  parameter-standard  deviation  of  material  concentration  m  the  Z 
direction  for  puff  of  interest 

Up  =  Wind  speed  of  individual  puff  at  the  time  of  release  (distance/time)*  at  10m 
tp  =  Time  since  puff  was  released  from  source 

H  =  Effective  stack  height  (includes  the  actual  stack  height  +  plume  rise  fi'om 
buoyant  forces) 

*  Time  units  for  source  term  (Q)  must  be  consistent  with  time  imits  for  wind  speed  (u) 
**  The  puff  is  assumed  to  have  a  circular  horizontal  cross  section;  therefore,  o  x  =■  ^  y 
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TABLE  4  Origin  of  Coordinates  for  Experiments  Conducted  During  Each  1  Hour  Period 


Treatment  1 

Treatment  2 

Treatment  3 

1-(1:1) 

10-(1;1) 

1-(1:2) 

10-(1:2) 

60-(l:2) 

1-(1:3) 

10-(1:3) 

60-(l:3) 

1-(10;1) 

10-(10:1) 

1-(10:2) 

10-(10:2) 

60-(10:2) 

1-(10:3) 

10-(10:3) 

60-(10;3) 

l-(60-l) 

10-(60:1) 

60-(60:l) 

l-(60-2) 

10-(60;2) 

l-(60-3) 

10-(60:3) 

“The  question  of  central  interest  here  is  whether  there  are  differences  in  true  averages 
associated  with  the  different  treatments . .  .’’(Devore,  1991)  In  this  experiment  the  null 
hypothesis  is  that  there  is  no  difference  between  plume  concentrations  in  the  receptor  grid 
when  one,  ten,  and  sixty  minute  weather  data  averaging  times  are  used.  The  alternate 
hypothesis  is  that  there  is  a  difference  between  plume  concentrations  in  the  receptor  grid 
when  one,  ten,  and  sixty  minute  weather  data  averaging  times  are  used.  If  the  ANOVA 
rejects  the  null  hypothesis,  Tukey’s  Procedure  is  used  to  determine  which  of  the 
population  means  are  different  from  one  another. 

STATISTDC  4.0  was  used  to  analyze  the  data  to  determine  whether  the  null 
hypothesis  was  true.  In  this  thesis  effort,  the  objective  was  to  study  whether  there  was  a 
statistically  significant  difference  in  plume  concentrations  modeled  using  one,  ten,  or  sixty 
minute  averaging  times  for  meteorological  data. 
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IV.  Results 


Chapter  Overview 

The  following  sections  examine  the  difference  between  plume  concentrations  when 
varying  averaging  times  are  used.  Analysis  of  Variance  (ANOVA)  and  Tukey’s  Procedure 
with  a  0.05  experiment- wise  error  rate  were  used  to  determine  whether  there  was  a 
statistically  significant  difference  between  the  use  of  one  minute,  ten  minute  or  sixty 
minute  meteorological  data  averaging  times  to  determine  plume  concentrations  in  a 
receptor  grid. 

Comparison  of  Maximum  Concentrations  by  Averaging  Time 

A  subjective  comparison  of  the  maximum  plume  concentrations  calculated  for  each 
of  the  three  averaging  times  was  accompHshed.  Table  5  represents  a  sample  of  the 
maximum  concentrations  achieved  in  the  receptor  ^d  dtiring  each  e)q)eriment.  For  each 
experimental  period  (1  hour),  the  maximum  plume  concentration  occurred  when  1  minute 
averaging  time  was  used.  The  60  minute  averaging  time  always  produced  the  smallest 
maximum  concentration.  A  probable  reason  for  this  result  is  the  calculation  used  to  derive 
the  00  values  for  10  minute  and  60  minute  averaging  times.  The  use  of  averaging  (use  of 
the  sample  mean)  to  compute  the  ten  minute  and  sixty  minute  weather  data  makes  the  data 
sensitive  to  outlying  values.  Confutation  of  the  value  of  (oq)  directly  effects  the  value  of 
the  horizontal  dispersion  coefficients  Ox  and  Oy.  Examination  of  the  Gaussian 
Instantaneous  Plume  Equation  (Eqn  1)  indicates  that  larger  values  of  Ox  and  Oy  yields 
smaller  plume  concentrations  Another  explanation  for  the  results  of  the  subjective 
maximum  concentration  comparison  is  the  speed  of  dispersal  of  the  plume.  The  1  minute 
averaging  time  ae  values  are  likely  to  be  smaller  than  the  60  minute  averaging  time  09  due 
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TABLE  5  Example  of  Maximum  Concentration  Reached  Within  the  Receptor  Grid  Based 

Upon  Specified  Averaging  Time  (Day  153) 


Time 

1  Minute 

10  Minute 

60  Minute 

1000 

1 

1 

5.26E-04 

3.02E-04 

L63E-04 

1100 

1 

5.14E-04 

2.83E-04 

L47E-04 

1200 

6.44E-04 

2.42E-04 

L44E-04 

1300 

7.29E-04 

2.81E-04 

L30E-04 

1400 

8.32E-04 

3.24E-04 

L67E-04 

1500 

3.93E-04 

L95E-04 

6.44E-05 

1600 

3.21E-04 

2.17E-04 

L12E-04 

1700 

3.77E-04 

2.35E-04 

1.28E-04 

1800 

7.57E-04 

2.71E-04 

1.49E.04 

to  the  averaging  used  to  derive  the  60  minute  ae.  Larger  values  of  a©  mean  larger  values 
of  Ox  and  Gy.  Larger  values  of  the  horizontal  dispersion  parameter  yield  smaller 
concentrations;  therefore,  the  1  minute  averaging  time  will  yield  larger  maximum 
concentrations. 

Comparison  of  Average  Concentrations  by  Averaging  Time 

A  subjective  comparison  of  the  average  concentrations  calculated  at  the  locations 
of  the  maximum  concentrations  show  that  the  average  concentration  using  sixty  minute 
averaging  time  were  larger  than  the  average  concentrations  at  one  minute  and  sixty  minute 
averaging  times.  In  aU  but  one  case  of  the  nine  examined,  the  sixty  minute  averaging  time 
yielded  the  highest  average  concentration  at  the  grid  point  of  the  largest  maxhnum 
concentration.  The  probable  cause  of  this  result  is  that  the  sixty  minute  meteorological 
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TABLE  6  Example  of  Average  Concentration  Reached  at  the  Location  of  the  Maximum 
Concentration  Within  the  Receptor  Grid  Based  Upon  Specified  Averaging  Time 


Time 

1 

! 

1  Minute 

10  Minute 

60  Minute 

1000 

L004E-04 

9.563E-05 

L589E-04 

1100 

2.383E-05 

1.030E-04 

1.431E-04 

1200 

6.383E-05 

7.487E-05 

4.104E-04 

1300 

4.083E-05 

8.700E-05 

1.258E-04 

1400 

4.863E-05 

1.161E-04 

1.617E-04 

1500 

2.607E-05 

5.847E-05 

6.317E-05 

1600 

4.825E-05 

9.032E-05 

L091E-04 

1700 

7.766E-05 

L042E-04 

L244E-04 

1800 

9.144E-05 

1.144E-04 

1.434E-04 

data  assumes  only  one  wind  speed  and  one  wind  direction  for  each  of  the  puffs  released; 
therefore,  they  foUow  one  another  on  their  trajectory  in  the  grid.  When  using  one  minute 
weather  data,  there  exists  the  possibility  that  sixty  separate  trajectories  are  generated  for 
each  puff  which  would  tend  to  spread  out  the  released  material  over  a  wider  area.  The 
same  is  true  often  minute  averaging,  but  release  would  likely  be  over  a  area  wider  than 
the  sixty  minute  data  area,  but  narrower  than  the  one  minute  data  area. 

Results  of  ANOVA  Based  on  Averaging  Times 

The  first  experiment  examined  using  ANOVA  and  Tukey’s  procedure  was  a  one- 
to-one  comparison  between  the  locations  of  the  three  maximum  concentrations  achieved 
in  the  receptor  grid  for  the  three  averaging  times  examined.  Table  6  indicates  the  total 
number  of  times  that  the  indicated  comparison  occurred.  This  table  indicates  that  the 
receptor  grid  concentrations  computed  using  one  minute  and  ten  minute  weather  data  are 
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TABLE  7  Total  Number  of  Statistically  Significant  Difierences  of  Calculated  Maxiinum 
Concentration  Between  the  Averaging  Times  Indicated 


Averaging  Times 

1  -  10 

10-60 

1-60 

Compared 

Differ 

Differ 

Differ 

Total 

15/26 

17/26 

26/26 

statistically  significantly  different  about  the  same  number  of  times  they  are  statistically  the 
same.  The  comparison  of  the  receptor  grid  concentrations  computed  using  ten  minute  and 
sixty  minute  weather  data  were  statistically  significantly  different  from  each  other 
approximately  65%  of  the  time  during  these  ejq)eriment.  However,  concentrations 
computed  using  one  minute  and  sixty  minute  weather  data  were  statistically  significantly 
different  in  all  cases  examined  by  this  research. 

Results  of  ANOVA  Based  on  Coordinates  Obtained  from  Averaging  Times 

The  second  experiment  conducted  using  ANOVA  and  Tukey’s  Procedure  was  a 
comparison  of  concentrations  at  coordinates  obtained  from  three  maximmn  concentrations 
for  each  averaging  time.  Table  7  illustrates  the  results  of  this  e>q)eriment.  For  exan[q)le, 
5/27  (row  2  -  column  2)  indicates  that  of  the  27  experiments  performed,  5  showed  a 
statistically  significant  difference  between  using  one  minute  and  ten  minute  weather  data  at 
the  coordinates  of  the  three  largest  concentrations  occurring  in  the  receptor  grid  when  one 
minute  averaging  time  was  used  to  determine  the  coordinates. 
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TABLE  8  Total  Number  of  Statistically  Significant  Differences  Between  the 
Concentrations  at  the  Coordinates  of  Averaging  Times  Indicated 


Origin  of  Coordinates 

j 

1-10 

1-60 

10-60 

(Averaging  Time) 

Differ 

Differ 

Differ 

1  Minute 

1 

5/27  i 

12/27 

4/27 

10  Minute 

13/27 

21/27 

9/27 

60  Minute 

18/26 

26/26 

17/26 

Using  coordinates  determined  from  the  locations  of  the  three  largest 
concentrations  using  one  minute  weather  data,  the  results  of  the  ANOVA  show  that  there 
in  these  experiments,  there  were  never  more  than  50%  of  the  experiments  showing  a 
statistically  significant  difference  between  averaging  times.  The  coordinates  from  the  10 
minute  weather  data  show  that  in  slightly  more  than  50%  of  the  experiments  show  a 
statistically  significant  difference  between  using  1  minute  and  10  minute  weather  data. 
However,  the  same  coordinates  show  that  in  more  than  75%  of  the  experiments,  there  was 
a  statistically  significant  difference  between  1  minute  and  60  minute  weather  data. 
Additionally,  the  coordinates  from  the  60  minute  meteorological  data  show  that  regardless 
of  averaging  time,  at  least  70%  of  the  experiments  show  a  statistically  significant 
difference  between  the  averaging  times.  Moreover,  the  60  minute  versus  1  minute 
weather  data  showed  a  statistically  significant  difference  in  100%  of  the  cases  examined. 

Summary 

In  each  of  the  experiments  performed  in  this  study,  100%  of  the  experiments 
showed  a  statistically  significant  difference  between  the  use  of  one  minute  versus  sixty 
minute  meteorological  data. 
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Comparison  of  ten  minute  averaging  time  to  one  minute  or  sixty  minute  averaging 
time  was  statistically  significantly  dififerent  in  more  than  50%  of  the  experiments. 
Although  when  comparisons  of  averaging  times  were  done  using  coordinates  generated 
fi:om  the  locations  of  the  three  largest  concentrations  for  each  averaging  time,  the  results 
show  that  less  than  50%  of  the  experiments  conqjaring  the  ten  minute  averaging  time  to 
the  one  minute  or  sixty  minute  averaging  time  are  statistically  significantly  different. 
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V.  Conclusions 


Introduction 

Analysis  of  the  receptor  grid  plume  concentrations  for  the  varying  averaging  times 
reveals  a  distinct  difference  between  using  one  minute  and  sixty  minute  averaging  times  in 
all  cases  examined  by  this  research  effort.  The  analysis  also  showed  differences  exist 
between  using  ten  minute  averaging  time  versus  one  minute  or  sixty  minute  averaging 
times.  This  chapter  translates  the  recognition  of  these  differences  among  averaging  times 
into  recommendations  for  future  modehng  of  Explosive  Ordnance  Disposal  (EOD)  Open 
Burning  (OB)  operations. 

Recommendations 

Maximum  Concentrations.  The  data  from  this  research  showed  that  in  all  cases, 
the  largest  maximum  concentrations  reached  in  each  receptor  grid  occurred  when  one 
minute  meteorological  data  was  used.  The  mostly  likely  cause  of  this  is  that  the  ten 
minute  and  sixty  minute  standard  deviations  of  the  wind  (sigma  thetas)  were  affected  by 
larger  outher  values  in  the  data  set.  A  remedy  for  this  problem  would  be  the  use  of  a 
“10%  trimmed  mean”  which  would  rid  the  data  set  of  extreme  values  (whether  large  or 
small).  The  use  of  the  trimmed  mean  might  bring  the  data  closer  to  the  true  mean. 

This  result  might  also  indicate  that  when  examhung  acute,  short-term,  exposure,  it 
is  more  important  to  use  one  minute  or  ten  minute  meteorological  data.  The  one-minute 
and,  to  a  lesser  extent,  ten  minute  data  are  less  likely  to  have  a  large  (sigma  theta)  during 
the  tunes  when  OB  is  hkely  to  take  place,  i.e.,  during  the  day.  This  smaller  (sigma  theta) 
signifies  the  likelihood  that  the  puff  wiU  not  difiuse  out  as  quickly  as  a  puff  with  a  larger 
(sigma  theta),  i.e.  sixty  minute  weather  data.  Due  to  the  large  number  of  chemicals  which 
could  be  released  during  an  OB  event,  constituents  with  acute  exposure  health  standards 
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should  be  examined  with  the  information  concerning  the  maximum  concentrations  reached 
with  one  minute  weather  data. 

Avera2e  Concentrations.  The  data  from  this  research  shows  that  in  all  but  one 
case  examined,  the  sixty  minute  averaging  time  yielded  the  highest  average  concentrations 
at  the  grid  points  where  maximum  concentrations  were  reached  in  the  receptor  grid.  This 
information  would  tend  to  be  useful  when  assessing  chronic,  long-term,  exposures  to 
vaiiousl  chemical  constituents.  The  use  of  one  set  of  meteorological  data  for  the  time 
period  of  interest  would  be  the  conservative  approach  because  the  puffs  released  during  an 
EOD  Open  Burning  operation  would  all  follow  the  same  trajectory;  therefore,  the  average 
concentration  would  tend  to  be  larger  than  puffs  released  with  differing  trajectories,  i.e, 
more  area  over  which  to  release  material. 

Comparison  of  Averagine  Times.  Analysis  of  Variance  (ANOVA)  and  Tukey’s 
Procedure  were  used  to  determine  whether  a  statistically  significant  difference  existed 
between  the  use  of  meteorological  data  with  varying  averaging  times.  In  each  case  where 
one  minute  was  compared  to  sbcty  minute  averaging  time,  a  statistically  significant 
difference  existed.  This  difference  was  most  likely  due  to  the  method  of  determing  the 
average  (sigma  theta)  for  the  averaging  time.  Outhers  in  the  data  set  can  have  a  profound 
effect  upon  the  concentrations  in  the  receptor  grid.  The  average  concentrations  for  the 
receptor  points  examined  in  this  research  effort  were  larger  for  the  sixty  minute  averaging. 
The  use  of  sixty  minute  averaging  time  weather  data  tends  to  support  the  current  air 
pollution  modeling  practice  of  one  set  of  weather  data  during  the  time  period  of  interest. 
Although,  the  weather  data  used  in  the  PUFFY  model  is  sixty  minute  weather  data,  ie., 
the  time  period  of  interest,  versus  the  ten  minute  weather  data  used  in  most  air  pollution 
models  through  out  the  tune  period  of  interest,  i.e.,  one  hour  in  this  case. 

Depending  upon  the  reason  why  the  PUFFY  model  is  being  used,  i.e.,  acute  vs 
chronic  exposures,  it  may  be  appropriate  to  use  one  averaging  time  for  acute  exposure  and 
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another  averaging  time  for  chronic  exposures.  In  any  case  it  is  important  for  future  EOD 
Open  Burning  models  to  account  for  this  dijBference  when  designing  the  model. 

Further  Research 

Several  opportunities  exist  for  future  research  on  the  issue  of  models  for  Explosive 
Ordnance  Disposal  OB  operations.  The  model  developed  for  this  research  eJBfort 
employed  several  simplifying  limitations.  The  effective  plume  height  was  assumed  in  the 
PUFFY  model,  while  a  more  comphcated  model  would  calculate  the  effective  plume 
height  based  upon  meteorological  data  as  well  as  data  about  the  OB  operation  (the 
source). 

The  PUFFY  model  used  an  imaginary  source  term  to  calculate  the  plume 
concentrations.  A  more  complex  model’s  use  of  emission  data  generated  by  Dugway 
Proving  Ground  experiments  to  develop  emissions  factor  could  be  used  to  develop  a  more 
reahstic  somrce  term.  The  PUFFY  model  used  a  static  value  of  the  source  term,  i.e.,  the 
same  source  strength  for  each  pufi^  a  more  complex  model’s  use  of  randomized  emission 
data  woixld  also  make  the  model  more  realistic.  The  time  of  release  for  the  puffs  was  set 
at  one  per  minute  for  the  PUFFY  model.  The  use  of  randomized  puff  release  times  would 
also  serve  to  make  the  model  mirror  reahty  more  closely. 

Finally,  this  research  effort  did  not  address  the  effects  of  the  time  of  day  on  the 
plume  concentrations  in  the  receptor  grid.  The  focus  of  fiiture  research  on  the  time  of  day 
would  seek  to  determine  whether  there  is  an  optimal  time  of  day  to  perform  an  EOD  Open 
Burning  operation.  This  could  possibly  help  impermitted  facUities  to  be  issued  a  permit  to 
bum  during  specified  hours  of  the  day. 
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Appendix  A.  ‘PUFFY’  FORTRAN  77  Code 


PROGRAM  PUFFY 
C 

C  FILE  NAME:  ' PUFFY. F' 

C 

C  DECLARATION  OF  VARIABLE  TYPES 
REAL  ALPHA,  STEP 

REAL  STRM,  XTRM,  YTRM,  ZTRM,  SUM,  CONCEN,  AGAMMA,  X,  Y,  HYP,  H 

INTEGER  TPER,  P,  TM,  PUFF,  REM,  XREF,  YREF,  GRIDX,  GRIDY 
INTEGER  HOUR,  HOURIO,  HOURS 0,  TIM,  MININCR,  INCR,  PERIOD,  TPERIOD 
INTEGER  XCOORD,  YCOORD,  N 
C 
C 

C  DIMENSION  STATEMENTS 

REAL  TEMP(60,4),  RH(60,4),  WSPEED (60 , 4) ,  WDIR(60,4) 

REAL  STHETA(60,4) 

REAL  DIST(60) ,TOP10 (10) ,COORDX(10) ,COORDY(10) ,  WHEN(IO) 

REAL  CONC(-101:101, -101:101) ,  MAXC ( - 101 : 101 , - 1 0 1 : 101 ) 

REAL  CONCMX(-101:101, -101:101) ,  AVG ( - 101 : 101 , - 101 : 101) 

REAL  SP(60),  A(60) ,  SOURCE(60) 

REAL  DTDZ(60),  ZSP(60),  C(60),  D(60) 

INTEGER  STN (60) ,  YR(60),  DATE(60),  TIME(60) 

INTEGER  HGT(60,4) 

C 

C 

C  DATA  STATEMENT 

Q* ************  *********************************************************  * 

C  SOURCE  -  Number  of  puffs  released  and  the  source  strength  assigned  to 

C  each  puff 

DATA  SOURCE/60*1000.0/ 

C 

C 

C  VARIABLE  INITIALIZATION 

REM  =  0 

TIM  =  0 

MININCR  =  0 
N  =  0 

CONCEN  =0.0 
AGAMMA  =0.0 
X  =  0.0 
Y  =  0.0 
HYP  =0.0 
STRM  =0.0 
XTRM  =  0.0 
YTRM  =0.0 
ZTRM  =0.0 
SUM  =0.0 
C 
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C 


C  CONSTANTS 

(]^****************'****'*-****i*r'**************************'********'**'***’****** 

C  H  -  Effective  Plume  Height  (meters) 

C  GRIDX  -  Number  of  grid  points  extending  from  the  source  in  the  +/- 

C  GRIDY  -  Number  of  grid  points  extending  from  the  source  in  the  +/- 

C  STEP  -  Distance  between  grid  points  (meters) 

C  TPER  -  Weather  Data  Averaging  Time 

C  INCR  -  Number  of  divisions  of  one  minute  used  to  for  calculations 

C  HOUR  -  Time  of  interest  to  begin  reading  1  minute  weather  data 

C  HODRIO  -  Time  of  interest  to  begin  reading  10  minute  weather  data 

C  HOUR60  -  Time  of  interest  to  begin  reading  60  minute  weather  data 

C  PERIOD  -  Length  of  time  of  interest  for  run  (hours) 

Q-k *********************************************************************  * 

H  =  25.0 

GRIDX  =20 

GRIDY  =20 

STEP  =50.0 

TPER  =  1 

INCR  =  4 

HOUR  =  1600 

HOURIO  =  HOUR  +10 

HOUR60  =  HOUR  +  100 

PERIOD  =  1 


FILE  DECLARATIONS 

WEATHER  DATA  FILES  GENERATED  BY  'metall.f': 

First  character,  i.e.,  'S',  =  station 

Second  character,  i.e.,  '1',  =  station  number 

Third  -  Fourth  characters,  i.e.,  '10',  =  averaging  time  (TPER) 
Remaining  characters,  i.e.,  'min',  =  minute 

OPEN  (UNIT  =  7,  FILE  =  's21min',  STATUS  =  'OLD') 

OPEN  (UNIT  =  8,  FILE  =  's210min',  STATUS  =  'OLD') 

OPEN  (UNIT  =  9,  FILE  =  's260rain',  STATUS  =  'OLD') 


'CONCl'  -  store  location  (x,y)  and  concentration  for  each  minute 

' PTMAXl '  -  store  location  (x,y)  and  value  of  maximum  concentration 

ever  reached  during  the  time  period  of  interest 
(for  each  grid  point) 

'MAXI'  -  store  location  (x,y)  and  value  of  cumulative  concentration 

C  reached  during  the  time  period  of  interest 

C  (for  each  grid  point) 

C  'AVGl'  -  store  location  (x,y)  and  value  of  average  concentration 

C  at  each  grid  point  during  the  time  period  of  interest 

Q*********************************************************************** 
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OPEN 

(UNIT 

= 

10, 

FILE 

= 

' concl '  , 

STATUS 

= 

'UNKNOWN' 

') 

OPEN 

(UNIT 

11, 

FILE 

= 

' ptmaxl ' , 

STATUS 

= 

'UNKNOWN' 

) 

OPEN 

(UNIT 

12, 

FILE 

= 

'  maxi '  , 

STATUS 

= 

'UNKNOWN' 

) 

OPEN 

(UNIT 

= 

13, 

FILE 

= 

'  avgl '  , 

STATUS 

= 

' UNKNOWN ' 

) 

OPEN 

(UNIT 

= 

14, 

FILE 

= 

' concl 0 ' , 

STATUS 

= 

' UNKNOWN ' 

) 

OPEN 

(UNIT 

= 

15, 

FILE 

= 

' ptmaxl 0 ' , 

STATUS 

= 

' UNKNOWN ' 

) 

OPEN 

(UNIT 

16, 

FILE 

= 

' maxlO ' , 

STATUS 

= 

'UNKNOWN' 

) 

OPEN 

(UNIT 

= 

17, 

FILE 

= 

' avglO ' , 

STATUS 

= 

' UNKNOWN ' 

) 

OPEN 

(UNIT 

= 

18, 

FILE 

' conc60 ' , 

STATUS 

= 

' UNKNOWN ' 

) 

OPEN 

(UNIT 

= 

19, 

FILE 

= 

' ptmax6  0 ' , 

STATUS 

= 

'UNKNOWN' 

) 

OPEN 

(UNIT 

20, 

FILE 

= 

' max6  0 ' , 

STATUS 

= 

' UNKNOWN ' 

) 

OPEN 

(UNIT 

= 

21, 

FILE 

' avg6  0 ' , 

STATUS 

= 

' UNKNOWN ' 

) 

OPEN 

(UNIT 

22, 

FILE 

= 

'  tenl '  , 

STATUS 

= 

' UNKNOWN ' 

) 

OPEN 

(UNIT 

= 

23, 

FILE 

= 

'tenlO' , 

STATUS 

= 

' UNKNOWN ' 

) 

OPEN 

(UNIT 

= 

24, 

FILE 

= 

' ten60 ' , 

STATUS 

' UNKNOWN ' 

) 

OPEN 

(UNIT 

= 

25, 

FILE 

= 

' statl ' , 

STATUS 

= 

' UNKNOWN ' 

) 

OPEN 

(UNIT 

= 

26, 

FILE 

= 

'statlO' , 

STATUS 

= 

' UNKNOWN ' 

) 

OPEN 

(UNIT 

= 

27, 

FILE 

' stateo '  , 

STATUS 

= 

' UNKNOWN ' 

) 

C 

C************************uSED  DURING  ANALYSIS  RUNS********************** 

C***************  take  comments  out  during  the  analysis  runs************* 
C  N  LOOP  -  used  to  read  the  appropriate  data  file  for  the  coordinates 
C  of  the  TOP  TEN  concentrations  reached  in  the  receptor  grid 

for  each  of  the  averaging  time 

DO  800  N  =  1,10 

IF  (TPER.EQ.l)  THEN 

READ  (22,*)  WHEN(N),  COORDX(N),  COORDY(N),  TOPIO(N) 

XCOORD  =  COORDX(N) 

YCOORD  =  COORDY(N) 

WRITE {*,*)  XCOORD , YCOORD 
ELSEIF  (TPER.EQ.IO)  THEN 

READ  (23,*)  WHEN(N),  COORDX(N),  COORDY{N),  TOPlO (N) 

C  XCOORD  =  COORDX(N) 

C  YCOORD  =  COORDY(N) 

C  WRITE(*,*)  XCOORD , YCOORD 

C  ELSEIF  (TPER.EQ.60)  THEN 

C  READ  (24,*)  WHEN(N),  COORDX(N),  COORDY(N),  TOPlO (N) 

C  XCOORD  =  COORDX(N) 

C  YCOORD  =  COORDY(N) 

C  WRITE (*,*)  XCOORD, YCOORD 

C  ENDIF 

C*************************  end  of  section*  *  ****************************** 

C  TPERIOD  -  time  interval  of  interest  (minutes) 

C  P  LOOP  -  Keeps  track  of  the  puff  releases  and  total  time  of  interest 
C  REM  -  used  to  determine  whether  proper  weather  data  is  read 

TPERIOD  =  PERIOD  *  60 
DO  600  P  =  1,  TPERIOD 
REM  =  MOD (P, 10) 
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c 

c 

Q*************************************************’********************** 

C  SUBROUTINE  initializes  the  concentration  grids,  i.e., 

C  the  changing  concentration  grid  and  the  grid  keeping  track  of  the 
C  maximum  concentration  at  each  grid  poing  during  a  2  hour  period 
(^•**'***************************************'*-*'************'************'**** 
CALL  INITIAL  (CONC,  GRIDX,  GRIDY,  P,  MAXC,  CONCMX) 

C 

IF  (N.GT.l)  GOTO  2000 
IF  (P.GT.60)  GOTO  1000 
C 

Q**********************************************-********’********-**'******'* 

C  READING  IN  THE  WEATHER  DATA  FOR  THE  TAGGED  PUFF 

C  ELSEIF  statements  are  used  to  make  sure  correct  meteorological  data 
C  is  read  in  at  the  correct  time,  i.e.,  one  puff  at  a  time 

C  TPER  =  1  -  new  data  will  be  read  in  each  time  through  the  loop 

C  TPER  =  10  -  new  data  will  be  read  in  after  each  set  of  10  runs 

C  TPER  =  60  -  same  data  for  each  puff 

C  RPTDATA  -  Reads  in  puff  data  that  is  the  same  as  previous  puff  data 

ie  ic  ic  i(  it  ic  it  ie  ie  ic  ie  -k  -k  it  ic  ic  -k  if  it  -k  ic  it  ie  it  "k  "k  ie  -k  ie  ic  ie  -k  ie  ie  "k  "ie  ic  ie  ie  ic  -k  ic  ir  ie  ie  "k  ie  ie  ic  it  ic  it  -k  ie  -k  -k  ic  it  "k  -k  ir  ie  ie  *  ie  ic  ie  ic  ic  "k 

100  IF  (TPER.EQ.l)  THEN 

READ(7,*)  STN(P),  YR(P),  DATE(P),  TIME(P), 

+  (HGT(P, J) ,TEMP (P, J) ,  RH(P,J),  WSPEED(P,J), 

+  WDIR{P,J),  STHETA(P, J) , J=l,4) 

C 

ELSEIF  (  (TPER.EQ.IO)  .AND.  (REM.EQ.l) )  THEN 
READ (8,*)  STN(P),  YR(P),  DATE{P),  TIME(P), 

+  (HGT(P, J) ,TEMP{P, J) ,  RH{P,J),  WSPEED(P,J), 

+  WDIR(P,J),  STHETA(P, J) , J=l,4) 

C 

ELSEIF  ( (TPER.EQ.IO) .AND. (REM. NE.l) )  THEN 
CALL  RPTDATA  (STN,  YR,  DATE,  TIME,  SP,  A,  DTDZ,  ZSP,  HGT, 

+  TEMP,  RH,  WSPEED,  WDIR,  STHETA,  P) 

GOTO  1000 
C 

ELSEIF  ( (TPER.EQ.60) .AND. (P.EQ.l) )  THEN 
READO,*)  STN(P),  YR(P),  DATE(P),  TIME(P), 

+  (HGT(P,J) ,TEMP(P, J) ,  RH(P,J),  WSPEED(P,J), 

+  WDIR(P,J),  STHETA(P, J) , J=l,4) 

C 

ELSEIF  ( (TPER.EQ.60) .AND. (P. NE.l) )  THEN 
CALL  RPTDATA  (STN,  YR,  DATE,  TIME,  SP,  A,  DTDZ,  ZSP,  HGT, 

+  TEMP,  RH,  WSPEED,  WDIR,  STHETA,  P) 

GOTO  1000 
END  IF 
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on  on 


([^*******************************************-jlr***'****-******************** 

C  These  statements  are  used  to  make  sure  the  data  is  read  from  the 
C  appropriate  place  in  the  data  file 

Qilf*'******’******************-*-******************************’************** 

IF  ( (TIME (P) .EQ.HOURIO) .AND. (P.EQ.l) )  GOTO  2000 
C 

IF  ( (TIME (P) .EQ.HOUR60) .AND. (TPER.EQ. 60) .AND. (P.EQ.l) )  THEN 
GOTO  2000 

ELSEIF  (P.EQ.l)  THEN 
GOTO  100 
END  IF 
C 

IF  ( (TIME (P) .EQ. HOUR) .AND. (TPER.NE. 10) .AND. (P.EQ.l) )  THEN 
GOTO  2000 

ELSEIF  (P.EQ.l)  THEN 
GOTO  100 
END  IF 
C 
C 

C  ********pRiNT  OUT  DATA  LINES  HERE  TO  CHECK  FOR  ACCURACY  ************** 
2000  WRITE(*,10)  STN(P),  YR(P),  DATE(P),  TIME(P), 

+  (HGT(P, J) ,TEMP(P, J) ,  RH(P,J),  WSPEED(P,J), 

+  WDIR(P,J),  STHETA(P, J) , J=l,4) 

10  F0RMAT(I2,I5,I4,I5,1X,4 (I2,3F5.1,  F6.1,F5.1)) 

WRITE (*,11)  P,  N 

11  FORMAT(2I4) 

^★★★★Vr***!!?***********************************'****************’************ 

CONVERSION  OF  WIND  SPEED  FROM  MPH  TO  M/S 

WSPEED(P,1)  =  WSPEED(P,l)*(1609./l.)*(l./60.)*(l./60.) 


C  SUBROUTINE  calculates  horizontal  and  vertical  stability  parameters 
C  and  horizontal  coefficients  and  exponents 

Q^^^lt********************-**********************************’*****’********* 

CALL  STABPAR(STHETA,  TEMP,  ZSP,  P,  A,  B) 

C 

C 

C******************************LOOP*********************************** 

C  TM  -  time  counter  helps  to  ensure  the  proper  puff  is  assigned 
C  the  proper  "age"  in  the  grid 

C  PUFF  -  tag  for  the  puff  of  interest 

Q************************Tkf******************************************** 

1000  DO  500  TM  =  1,P 

PUFF  =  P  -  TM  +  1 
IF  (PUFF.GT.60)  GOTO  500 
C 
C 
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C*  ************************  *****l^QQp***  ******************************  ***■), 

C  MININCR  -  this  counter  is  set  to  the  number  of  partitions  desired 
C  per  minute,  i.e.,  (1,4  =>  15  minute  increments), 

C  (1,10  =>  6  minute  increments) 

C  (Depending  upon  what  MININCR  is  changed  to,  the  variable  TIM  must 
C  also  be  changed  (automatically)  to  reflect  the  new  increment,  i.e., 

C  [15  minute  increments  =>  TIM  =  (TM*60)  +  (MININCR*15)  -  60] , 

C  [6  minute  increments  =>  TIM  =  (TM*60)  +  (MININCR*6)  -  60] ) 

C  DIST  -  distance  downwind  from  source  to  center  of  puff 
C  (SPEED*TIME  =  DISTANCE) 

DO  400  MININCR  =  1 , INCR 

TIM  =  (TM*60)  +  (MININCR* (60/INCR) )  -  60 
DIST (PUFF)  =  WSPEED (PUFF, 1) * (TIM) 

C 

************************************************************************* 

SUBROUTINE  calculates  vertical  diffusion  coefficients  and  exponents 
*********•***★*■*■******■*★**■*■****•******★★*★★***★★************■************* 

CALL  CDSTAB(2SP,  DIST,  C,  D,  PUFF) 


SUBROUTINE  calculates  wind  direction  in  +  Cartesian  Coordinates 
CALL  DEGREE  (WDIR,  ALPHA,  P) 


★  ★★★★★★★★★★★★★★★★★★★★★if*************************’********************** 

SUBROUTINE  calculates  the  concentrations  at  each  point  in  the 
stationary  frame  of  reference  grid  at  each  individual  minute 

CALL  REFGRID(XREF,  YREF,  GRIDX,  GRIDY,  ALPHA,  PUFF,  A,  B,  C,  D, 

+  DIST,  SOURCE,  WSPEED,  TIM,  CONC,  STEP,  H,  INCR) 

GO  BACK  TO  CHANGE  THE  MINUTE  INCREMENT  (MININCR) 

400  CONTINUE 

C 

C  GO  BACK  TO  CHANGE  THE  TIME  INCREMENT  (TM) 

500  CONTINUE 
C 
C 

C  SUBROUTINE  calculates  the  maximum  concentration  ever  reached  during 
C  the  time  period  of  interest  (at  each  grid  point) 

CALL  PTMAX  (CONC,  CONCMX,  GRIDX,  GRIDY) 

C 

C 
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C  SUBROUTINE  calculates  the  cumulative  concentration  at  each  point  in 
C  the  stationary  frame  of  reference  grid 

CALL  MAXCONC  (CONC,  MAXC,  GRIDX,  GRIDY) 

^ic  ic  ic  ic  ie  ic  ie  it  ic  ic  ic  ic  -k  "k  "k  it  "k  -k  it  -k  it  ie  ic  "k  ic  if  ie  ic  "k  "k  "k  "k  -k  it  "k  ic  "k  ic  "k  ie  ic  ie  -k  -k  ie  "k  ie  ic  -k  ic  ic  ie  "k  "k  ic  -k  ie  ic  ie  -k  i(  ie  -k  ie  "k  -k  -k  ic  ic 

C  SUBROUTINE  calculates  the  "Top  Ten"  concentrations  reached  in  the 
C  receptor  grid 

Qk  k  k  k  it  k  it  it  it  it  k  k  k  it  k  k  it  k  k  k  k  k  k  k  k  k  it  it  k  k  k  k  k  k  k  k  k  k  k  it  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  it  k  k  it  k 

CALL  TOPTEN  (P,  GRIDX,  GRIDY,  CONC,  TOPIO,  COORDX,  COORDY,  WHEN) 
C 
C 

C************************uSED  DURING  REGULAR  RUNS********************* 
C*****C0MMENT  (C)  THIS  SECTION  OUT  WHEN  DOING  THE  ANALYSIS  RUN******** 
C  Writing  the  reference  grid  to  the  screen  before  returning  to  read 
C  data  about  another  puff 

C  The  reference  grid  is  sized  to  fit  the  screen 

C  Writes  locations  and  concentrations  to  a  file  for  later  graphing 
C  File  Name:  'concl' 

C 

IF  (TPER.EQ.l)  THEN 

WRITE (*,*)  'CONCENTRATION  AT  EACH  POINT  DURING  THIS  RUN' 
WRITE(*,19)  ( (CONC{XREF,YREF) ,  XREF  =  -6,6), 

+  YREF=GRIDY, -GRIDY, -1) 

19  FORMAT (13E10. 3) 

DO  110  XREF  =  -GRIDX, GRIDX, 1 
DO  120  YREF  =  GRIDY, -GRIDY, - 1 

WRITE(10,23)  XREF,  YREF,  CONC (XREF, YREF) 

120  CONTINUE 

110  CONTINUE 
23  FORMAT(2I4,  E10.3) 

WRITE (*,21)  P 
WRITE (10,21)  P 
21  FORMAT (15) 

C 

ELSEIF  (TPER.EQ.IO)  THEN 

WRITE (*,*)  'CONCENTRATION  AT  EACH  POINT  DURING  THIS  RUN' 

WRITE (*,19)  ( (CONC (XREF, YREF) ,  XREF  =  -6,6), 

+  YREF=GRIDY, -GRIDY, -1) 

DO  130  XREF  =  -GRIDX, GRIDX 

DO  140  YREF  =  GRIDY, -GRIDY, - 1 

WRITE(14,23)  XREF,  YREF,  CONC (XREF, YREF) 

140  CONTINUE 

130  CONTINUE 

WRITE (*,21)  P 
WRITE (14,21)  P 
C 

ELSEIF  (TPER.EQ.60)  THEN 

WRITE (*,*)  'CONCENTRATION  AT  EACH  POINT  DURING  THIS  RUN' 
WRITE(*,19)  ( (CONC (XREF, YREF) ,  XREF  =  -6,6), 

+  YREF=GRIDY, -GRIDY, -1) 


nnnnnnnnnn 


DO  150  XREF  =  -GRIDX,GRIDX 

DO  160  YREF  =  GRIDY, -GRIDY, -1 

WRITE(18,23)  XREF,  YREF,  CONG (XREF, YREF) 
160  CONTINUE 

150  CONTINUE 

WRITE (*,21)  P 
WRITE(18,21)  P 
END  IF 


C;*************************enD  OF  SECTION******************************** 

C 

C 

C************************usED  DURING  ANALYSIS  RUNS********************** 
C**************taKE  comments  OUT  DURING  THE  ANALYSIS  RUNS*************** 
C 

C  IF  (TPER.EQ.l)  THEN 

C  WRITE (*,*)  'CONCENTRATION  AT  THIS  POINT  DURING  THIS  RUN' 

C  WRITE (*,20)  CONC(XCOORD,YCOORD) 

C  WRITE (25,20)  CONC (XCOORD, YCOORD) 

C  20  FORMAT (ElO. 3) 

C  WRITE (*,22)  P 

C  WRITE (25, 22)  P 

C  22  FORMAT (13) 

C  ELSEIF  (TPER.EQ.IO)  THEN 

C  WRITE (*,*)  'CONCENTRATION  AT  THIS  POINT  DURING  THIS  RUN' 

C  WRITE (*,20)  CONC (XCOORD, YCOORD) 

C  WRITE (26,20)  CONC (XCOORD, YCOORD) 

C  WRITE  (*,22)  P 

C  WRITE (26,22)  P 

ELSEIF  (TPER.EQ.60)  THEN 

WRITE (*,*)  'CONCENTRATION  AT  THIS  POINT  DURING  THIS  RUN' 

WRITE (*,20)  CONC (XCOORD, YCOORD) 

WRITE(27,20)  CONC (XCOORD, YCOORD) 

WRITE (*,22)  P 
WRITE  (27,22)  P 
END  IF 


******************** *****end  OF  SECTION******************************** 


C 

C  GO  BACK  TO  CALCULATE  INFORMATION  FOR  ANOTHER  PUFF 
600  CONTINUE 
C 
C 

C  SUBROUTINE  calculates  the  average  concentration  at  each  grid 
C  point  over  the  time  period  of  interest 

CALL  AVRG(AVG,  MAXC,  TPERIOD,  GRIDX,  GRIDY) 

C 

C 
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(;;;************************USED  DURING  REGULAR  RUNS*********************** 
C*****COMMENT  (C)  THIS  SECTION  OUT  WHEN  DOING  THE  ANALYSIS  RUN********* 

C  Writing  the  maximum  concentration  to  ever  occur  at  each  grid  point 
C  to  the  screen 
C  File  Name:  'ptmaxl' 

C 

IF  (TPER.EQ.l)  THEN 

WRITE  (*,*)  'MAXIMUM  CONCENTRATION  AT  EACH  POINT  AT  ANY  TIME' 
WRITE (*,19)  ( (CONCMX(XREF,YREF) ,  XREF  =  -6,6), 

+  YREF=GRIDY, -GRIDY, -1) 

DO  170  XREF  =  -GRIDX,GRIDX 

DO  180  YREF  =  GRIDY, -GRIDY, -1 

WRITE (11,23)  XREF,  YREF,  CONCMX (XREF, YREF) 

180  CONTINUE 

170  CONTINUE 

WRITE (*,*)  'TPER  =  ',  TPER 
WRITE (11,21)  TPER 
C 

ELSEIF  (TPER. EQ. 10)  THEN 

WRITE (*,*)  'MAXIMUM  CONCENTRATION  AT  EACH  POINT  AT  ANY  TIME' 
WRITE  (*,19)  ( (CONCMX (XREF, YREF)  ,  XREF  =  -6,6), 

+  YREF=GRIDY, -GRIDY, -1) 

DO  190  XREF  =  -GRIDX,GRIDX 

DO  210  YREF  =  GRIDY, -GRIDY, -1 

WRITE (15,23)  XREF,  YREF,  CONCMX (XREF, YREF) 

210  CONTINUE 

190  CONTINUE 

WRITE (*,*)  'TPER  =  ',  TPER 
WRITE(15,21)  TPER 
C 

ELSEIF  (TPER. EQ. 60)  THEN 

WRITE  (*,*)  'MAXIMUM  CONCENTRATION  AT  EACH  POINT  AT  ANY  TIME' 
WRITE(*,19)  ( (CONCMX (XREF, YREF) ,  XREF  =  -6,6), 

+  YREF=GRIDY, -GRIDY, -1) 

DO  220  XREF  =  -GRIDX,GRIDX 

DO  230  YREF  =  GRIDY, -GRIDY, -1 

WRITE (19,23)  XREF,  YREF,  CONCMX (XREF, YREF) 

230  CONTINUE 

220  CONTINUE 

WRITE  (*,*)  'TPER  =  ',  TPER 
WRITE (19,21)  TPER 
END  IF 


Q*’****'***'*'’************-}!?*'**'*’* END  OF  SECTION* ***********'*^'f**********'*r***** 

c 

c 
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C************************USED  DURING  ANALYSIS  RUNS********************** 
C**************i’AKE  COMMENTS  OUT  DURING  THE  ANALYSIS  RUNS*************** 
C  IF  (TPER.EQ.l)  THEN 

C  WRITE (*,*)  'MAXIMUM  CONCENTRATION  AT  THIS  POINT  AT  ANY  TIME' 

C  WRITE (*,20)  CONCMX(XCOORD,YCOORD) 

C  WRITE (*,*)  'TPER  =  ',  TPER 

C 

C  ELSEIF  (TPER. EQ. 10)  THEN 

C  WRITE (*,*)  'MAXIMUM  CONCENTRATION  AT  THIS  POINT  AT  ANY  TIME' 

C  WRITE (*,20)  CONCMX(XCOORD,YCOORD) 

C  WRITE(*,*)  'TPER  =  ',  TPER 

C 

C  ELSEIF  (TPER. EQ. 60)  THEN 

C  WRITE (*,*)  'MAXIMUM  CONCENTRATION  AT  THIS  POINT  AT  ANY  TIME' 

C  WRITE (*,20)  CONCMX (XCOORD, YCOORD) 

C  WRITE (*,*)  'TPER  =  ',  TPER 

C  ENDIF 

C 

C****************************enD  OF  SECTION***************************** 

C 

C 

C************************USED  DURING  REGULAR  RUNS*********************** 
C*****C0MMENT  (C)  THIS  SECTION  OUT  WHEN  DOING  THE  ANALYSIS  RUN********** 
C  Writing  the  cumulative  concentration  at  each  point  in  the  receptor 
C  grid  to  the  screen  and  to  a  file 
C  File  Name:  'maxi' 

C 

IF  (TPER.EQ.l)  THEN 

WRITE (*,*)  'CUMULATIVE  CONCENTRATION' 

WRITE (*,19)  ( (MAXC(XREF,YREF) ,  XREF  =  -6,6), 

+  YREF=GRIDY, -GRIDY, -1) 

DO  240  XREF  =  -GRIDX,GRIDX 

DO  250  YREF  =  GRIDY, -GRIDY, -1 

WRITE (12,23)  XREF,  YREF,  MAXC (XREF, YREF) 

250  CONTINUE 

240  CONTINUE 

WRITE (*,*)  'TPER  =  ',  TPER 
WRITE (12,21)  TPER 
C 

ELSEIF  (TPER. EQ. 10)  THEN 

WRITE (*,*)  'CUMULATIVE  CONCENTRATION' 

WRITE(*,19)  ( (MAXC (XREF, YREF) ,  XREF  =  -6,6), 

+  YREF=GRIDY, -GRIDY, -1) 

DO  260  XREF  =  -GRIDX,GRIDX 

DO  270  YREF  =  GRIDY, -GRIDY, -1 

WRITE(16,23)  XREF,  YREF,  MAXC (XREF , YREF) 

270  CONTINUE 

260  CONTINUE 

WRITE (*,*)  'TPER  =  ',  TPER 
WRITE(16,21)  TPER 
C 
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ELSEIF  {TPER.EQ.60)  THEN 

WRITE (*,*)  'CUMULATIVE  CONCENTRATION' 

WRITE(*,19)  ( (MAXC(XREF,YREF) ,  XREF  =  -6,6), 

+  YREF=GRIDY, -GRIDY, -1) 

DO  280  XREF  =  -GRIDX,GRIDX 

DO  290  YREF  =  GRIDY, -GRIDY, -1 

WRITE(20,23)  XREF,  YREF,  MAXC (XREF , YREF) 

290  CONTINUE 

280  CONTINUE 

WRITE (*,*)  'TPER  =  ’,  TPER 
WRITE (20,21)  TPER 
END  IF 
C 

0***************************gND  OF  SECTION****************************** 

C 

C 

c************************used  during  analysis  runs********************** 
C**********takE  out  comments  (C)  when  DOING  THE  ANALYSIS  RUN************ 
C 

C  IF  (TPER.EQ.l)  THEN 

C  WRITE (*,*)  'CUMULATIVE  CONCENTRATION' 

C  WRITE (*,20)  MAXC(XCOORD,YCOORD) 

C  WRITE (*,*)  'TPER  =  ',  TPER 

C 

C  ELSEIF  (TPER. EQ. 10)  THEN 

C  WRITE (*,*)  'CUMULATIVE  CONCENTRATION' 

C  WRITE (*,20)  MAXC(XCOORD,YCOORD) 

C  WRITE (*,*)  'TPER  =  ',  TPER 

C 

ELSEIF  (TPER.EQ.60)  THEN 

WRITE (*,*)  'CUMULATIVE  CONCENTRATION' 

C  WRITE (*,20)  MAXC(XCOORD,YCOORD) 

C  WRITE (*,*)  'TPER  =  ',  TPER 

C  ENDIF 

C 

c***************************end  of  section****************************** 

C 

C 

0************************USED  DURING  REGULAR  RUNS*********************** 
C*****COMMENT  (C)  THIS  SECTION  OUT  WHEN  DOING  THE  ANALYSIS  RUN********** 
C  Writing  the  average  concentration  at  each  point  in  the  receptor 
C  grid  to  the  screen  and  to  a  file 
C 

IF  (TPER.EQ.l)  THEN 

WRITE ( * , * )  ' AVERAGE  CONCENTRATION ' 

WRITE (*,19)  ( (AVG (XREF, YREF) ,  XREF  =  -6,6), 

+  YREF=GRIDY, -GRIDY, -1) 

DO  310  XREF  =  -GRIDX,GRIDX 

DO  320  YREF  =  GRIDY, -GRIDY, -1 

WRITE(13,23)  XREF,  YREF,  AVG (XREF, YREF) 

320  CONTINUE 

310  CONTINUE 
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WRITE (*,*)  'TPER  =  TPER 
WRITE (13,21)  TPER 
C 

ELSEIF  (TPER. EQ. 10)  THEN 

WRITE ( * , * )  ' AVERAGE  CONCENTRATION ' 

WRITE (*,19)  ( (AVG(XREF,YREF) ,  XREF  =  -6,6), 

+  YREF=GRIDY. -GRIDY, -1) 

DO  330  XREF  =  -GRIDX,GRIDX 

DO  340  YREF  =  GRIDY, -GRIDY, -1 

WRITE(17,23)  XREF,  YREF,  AVG (XREF, YREF) 

340  CONTINUE 

330  CONTINUE 

WRITE (*,*)  'TPER  =  TPER 
WRITE (17,21)  TPER 
C 

ELSEIF  (TPER.EQ.60)  THEN 

WRITE ( * , * )  ' AVERAGE  CONCENTRATION ' 

WRITE (*,19)  ( (AVG (XREF, YREF) ,  XREF  =  -6,6), 

+  YREF=GRIDY, -GRIDY, -1) 

DO  350  XREF  =  -GRIDX,GRIDX 
DO  360  YREF  =  GRIDY, -GRIDY, - 1 

WRITE (21, 23)  XREF,  YREF,  AVG (XREF, YREF) 

360  CONTINUE 

350  CONTINUE 

WRITE(*,*)  'TPER  =  ',  TPER 
WRITE (21,21)  TPER 
END  IF 
C 

C***************************enD  of  SECTION****************************** 

C 

C 

C************************USED  DURING  ANALYSIS  RUNS********************** 
C**********TAKE  OUT  COMMENTS  (C)  WHEN  DOING  THE  ANALYSIS  run************ 
C 

C  IF  (TPER.EQ.l)  THEN 

C  WRITE (*,*)  'AVERAGE  CONCENTRATION' 

C  WRITE (*,20)  AVG (XCOORD , YCOORD) 

C  WRITE (*,*)  'TPER  =  ',  TPER 

ELSEIF  (TPER. EQ. 10)  THEN 
C  WRITE (*,*)  'AVERAGE  CONCENTRATION' 

C  WRITE (*,20)  AVG (XCOORD, YCOORD) 

C  WRITE (*,*)  'TPER  =  ',  TPER 

C 

C  ELSEIF  (TPER.EQ.60)  THEN 

C  WRITE (*,*)  'AVERAGE  CONCENTRATION' 

C  WRITE (*,20)  AVG (XCOORD, YCOORD) 

C  WRITE (*,*)  'TPER  =  ',  TPER 

C  ENDIF 

C 

C**************************enD  of  SECTION******************************* 

C 
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C************************uSED  DURING  REGULAR  RUNS*********************** 
C*******C0MMENT  (C)  OUT  THIS  ENTIRE  SECTION  DURING  THE  ANALYSIS  RUN***** 
C  Writing  the  TOP  TEN  concentrations  to  the  screen  and  to  a  file 
C 

IF  (TPER.EQ.l)  THEN 

WRITE (*,*)  'TOP  TEN  CONCENTRATIONS  AND  WHEN/WHERE  THEY  OCCUR' 

DO  725  I  =  1,10 

WRITE (*,25)  WHEN(I),  COORDX(I),  COORDY(I),  TOPIO(I) 
WRITE(22,25)  WHEN(I),  COORDX(I),  COORDY(I),  TOPIO(I) 

25  FORMAT(3I5,E12 .3) 

725  CONTINUE 

WRITE(*,*)  'TPER  =  ',  TPER 
WRITE (22,21)  TPER 
WRITE (22,21)  HOUR 
C 

ELSEIF  (TPER. EQ. 10)  THEN 

WRITE (*,*)  'TOP  TEN  CONCENTRATIONS  AND  WHEN/WHERE  THEY  OCCUR' 

DO  750  I  =  1,10 

WRITE(*,25)  WHEN(I),  COORDX(I),  COORDY(I),  TOPIO(I) 

WRITE(23,25)  WHEN(I),  COORDX(I),  COORDY(I),  TOPIO(I) 

750  CONTINUE 

WRITE (*,*)  'TPER  =  ',  TPER 
WRITE (23,21)  TPER 
WRITE (23,21)  HOURIO 
C 

ELSEIF  (TPER. EQ. 60)  THEN 

WRITE ( * , * )  ' TOP  TEN  CONCENTRATIONS  AND  WHEN/WHERE  THEY  OCCUR ' 
WRITE ( * , * )  ' TOP  TEN  CONCENTRATIONS  AND  WHEN/WHERE  THEY  OCCUR ' 

DO  775  I  =  1,10 

WRITE(*,25)  WHEN(I),  COORDX(I),  COORDY(I),  TOPIO(I) 

WRITE(24,25)  WHEN(I),  COORDX(I),  COORDY(I),  TOPIO(I) 

775  CONTINUE 

WRITE (*,*)  'TPER  =  ',  TPER 
WRITE (24,21)  TPER 
WRITE (24,21)  HOUR60 
END  IF 
C 

C****** ************** ********EISfD  OF  SECTION***************************** 

C 

C 

(;;;************************USED  DURING  ANALYSIS  RUNS********************** 
C************coMMENT  THIS  STATEMENT  OUT  WHEN  DOING  REGULAR  RUN********** 
C  End  of  N  loop 
C 

C  800  CONTINUE 

(;;;***************************END  OF  SECTION****************************** 

C 

C 

c**************************end  of  main  program************************** 

STOP 

END 
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SUBROUTINE  AVRG(AVG,  MAXC,  TPERIOD,  GRIDX,  GRIDY) 

C 

C  THIS  SUBROUTINE  CALCULATES  THE  AVERAGE  CONCENTRATION  OVER  THE  TIME 
C  PERIOD  OF  INTEREST 
C 

C  DECLARATION  OF  VARIABLE  TYPES 

INTEGER  GRIDX,  GRIDY,  TPERIOD 
C 

C  DIMENSION  STATEMENTS 

REAL  AVG(-101:101,  -101:101),  MAXC ( -101 : 101 ,  -101:101) 

C 

DO  100  XREF  =  -GRIDX,  GRIDX 

DO  200  YREF  =  GRIDY,  -GRIDY,  -1 

AVG{XREF,YREF)  =  MAXC (XREF, YREF) /TPERIOD 
200  CONTINUE 
100  CONTINUE 
C 

C  WRITE  (*,15)  ( (AVG (XREF, YREF)  , XREF  =  -GRIDX, GRIDX)  , 

C  +  YREF  =  GRIDY, -GRIDY, -1) 

15  FORMAT  (13E10.3) 

RETURN 
END 
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SUBROUTINE  CDSTAB(ZSP,  DIST,  C,  D,  PUFF) 

C 

C  THIS  SUBROUTINE  CALCULATES  VERTICAL  STABILITY  PARAMETER  COEFFICIENTS 
C  AND  EXPONENTS 
C 

C  DECLARATION  OF  VARIABLE  TYPES 

INTEGER  PUFF 
C 

C  DIMENSION  STATEMENTS 

REAL  DIST(60) ,  C{60),  D(60),  ZSP(60) 

C 

IF  ( (ZSP(PUFF) .LE.1.0) .AND. (DIST(PUFF) .LT.745.) )  THEN 
C(PUFF)  =  0.0414 
D(PUFF)  =  1.3155 

ELSEIF  { (ZSP(PUFF) .LE.1.0) .AND. {DIST(PUFF) .GE. 745.) )  THEN 
C(PUFF)  =  1.928E-4 
D{PUFF)  =  2.1234 

ELSEIF  (  (ZSP(PUFF)  .LE. 2.0)  .AND.  (DIST(PUFF)  .LT.745.) )  THEN 
C(PUFF)  =  0.1036 
D(PUFF)  =  1.0026 

ELSEIF  ( (ZSP(PUFF) .LE.2.0) .AND. (DIST(PUFF) .GE.745.) )  THEN 
C(PUFF)  =  0.0534 
D{PUFF)  =  1.1029 

ELSEIF  {  (ZSP(PUFF)  .LE. 3.0)  .AND.  (DIST(PUFF)  .LT. 2000.) )  THEN 
C(PUFF)  =  0.1173 
D(PUFF)  =  0.9112 

ELSEIF  {  (ZSP(PUFF)  .LE. 3.0)  .AND. (DIST(PUFF)  .GE. 2000.) )  THEN 
C(PUFF)  =  0.4422 
D(PUFF)  =  0.7382 

ELSEIF  { (ZSP(PUFF) .LE. 4.0) .AND. (DIST(PUFF) .LT. 1100.) )  THEN 
C(PUFF)  =  0.0975 
D{PUFF)  =  0.8414 

ELSEIF  ( (ZSP(PUFF) .LE. 4.0) .AND. (DIST(PUFF) .GE. 1100.) )  THEN 
C(PUFF)  =  0.6097 
D(PUFF)  =  0.5808 

ELSEIF  ( (ZSP(PUFF) .LE. 5.0) .AND. (DIST{PUFF) .LT. 1400.) )  THEN 
C(PUFF)  =  0.1050 
D(PUFF)  =  0.7692 

ELSEIF  ( (ZSP(PUFF) .LE. 5.0) .AND. (DIST(PUFF) .GE. 1400.) )  THEN 
C(PUFF)  =  0.8788 
D(PUFF)  =  0.4771 

ELSEIF  ( (ZSP(PUFF) .GT. 5.0) .AND. (DIST(PUFF) .LT. 1400.) )  THEN 
C(PUFF)  =  0.0617 
D(PUFF)  =  0.7884 

ELSEIF  ( (ZSP(PUFF) .GT. 5.0) .AND. (DIST(PUFF) .GE. 1400.) )  THEN 
C(PUFF)  =  0.9990 
D(PUFF)  =  0.4771 

END  IF 

RETURN 

END 
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SUBROUTINE  DEGREE  (WDIR,  ALPHA,  P) 


C 

C  THIS  SUBROUTINE  CONVERTS  WIND  DIRECTION  TO  +  CARTESIAN  (DEGREES) 
C 

C  DECLARATION  OF  VARIABLE  TYPES 

REAL  DIRNEW,  XNEW,  YNEW,  XTRANS,  YTRANS,  ALPHA 
INTEGER  P 
C 

C  DIMENSION  STATEMENTS 
REAL  WDIR  (60,4) 

C 

DIRNEW  =  WDIR(P,1)  +  180. 

IF  (DIRNEW.lt. 0.0)  DIRNEW  =  DIRNEW  +  360. 

XNEW  =  COSD (DIRNEW) 

YNEW  =  SIND (DIRNEW) 

XTRANS  =  YNEW 
YTRANS  =  XNEW 

ALPHA  =  ATAN2D( YTRANS, XTRANS) 

IF  (ALPHA.lt. 0.0)  ALPHA  =  ALPHA  +  360. 

RETURN 

END 
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SUBROUTINE  INITIAL  (CONC,  GRIDX,  GRIDY,  P,  MAXC,  CONCMX) 
C  THIS  SUBROUTINE  ZEROES  OUT  THE  GRID  WHEN  NECESSARY 
C 

C  DECLARATION  OF  VARIABLE  TYPES 
INTEGER  GRIDX,  GRIDY,  P 
C 

C  DIMENSION  STATEMENTS 

REAL  CONC{-101;101, -101:101) ,  MAXC (- 101 : 101 ,- 101 : 101 ) 
REAL  CONCMX(-101:101, -101:101) ,  TOPIO(IO) 

INTEGER  COORDX( 10) ,  COORDY(IO) 

C 

DO  100  I  =  -GRIDX,  GRIDX 
DO  200  J  =  GRIDY, -GRIDY, -1 
CONC(I,J)  =0.0 
200  CONTINUE 
100  CONTINUE 
C 

IF  (P.GT.l)  GOTO  1100 
C 

DO  300  I  =  -GRIDX,  GRIDX 

DO  400  J  =  GRIDY,  -GRIDY,  -1 
CONCMX { I, J)  =0.0 
MAXC (I, J)  =0.0 
400  CONTINUE 
300  CONTINUE 
C 

DO  700  I  =  1,10 
TOPIO(I)  =  0.0 
COORDX(I)  =  0 
COORDY(I)  =  0 
700  CONTINUE 
C 

1100  RETURN 
END 
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SUBROUTINE  MAXCONC ( CONG ,  MAXC,  GRIDX,  GRIDY) 

C 

C  THIS  SUBROUTINE  CALCULATES  THE  MAXIMUM  CONCENTRATION  AT  EACH  GRID 
C  POINT  DURING  A  TWO  HOUR  PERIOD 
C 

C  DECLARATION  OF  VARIABLE  TYPES 

INTEGER  XREF,  YREF,  GRIDX,  GRIDY 
C 

C  DIMENSION  STATEMENTS 

REAL  CONC{-101:101,  -101:101),  MAXC ( -101 : 101 ,  -101:101) 

C 

DO  100  XREF  =  -GRIDX,  GRIDX 

DO  200  YREF  =  GRIDY,  -GRIDY,  -1 

MAXC (XREF, YREF)  =  MAXC (XREF, YREF)  +  CONC (XREF, YREF) 

200  CONTINUE 
100  CONTINUE 
C 

RETURN 

END 
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SUBROUTINE  PTMAX ( CONG ,  CONCMX,  GRIDX,  GRIDY) 

C 

C  THIS  SUBROUTINE  FINDS  THE  MAXIMUM  CONCENTRATION  AT  ANY  ONE  GRID  POINT 
C  DURING  THE  TWO  HOUR  PERIOD  OF  INTEREST 
C 

C  DECLARATION  OF  VARIABLE  TYPES 

INTEGER  XREF,  YREF,  GRIDX,  GRIDY 
C 

C  DIMENSION  STATEMENTS 

REAL  CONC(-101:101, -101:101) ,  CONCMX ( -101 : 101 ,  -101:101) 

C 

DO  100  XREF  =  -GRIDX,  GRIDX 

DO  200  YREF  =  GRIDY,  -GRIDY,  -1 
IF  (P.EQ.l)  GOTO  1000 
GOTO  1100 

1000  CONCMX (XREF, YREF)  =  CONC (XREF, YREF) 

1100  IF  ( CONC (XREF, YREF) .GT. CONCMX (XREF, YREF) ) 

+  CONCMX (XREF, YREF)  =  CONC (XREF, YREF) 

200  CONTINUE 
100  CONTINUE 
RETURN 
END 
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SUBROUTINE  REFGRID (XREF,  YREF,  GRIDX,  GRIDY,  ALPHA,  PUFF,  A,  B, 

+  C,  D,  DIST,  SOURCE,  WSPEED,  TIM,  CONC,  STEP,  H,  INCR) 

C 

C  DECLARATION  OF  VARIABLE  TYPES 

REAL  X,  Y,  Z,  SQ,  HYP,  BETA,  GAMMA,  AGAMMA,  ALPHA,  XSHIFT,  YSHIFT 
REAL  STRM,  XTRM,  YTRM,  ZTRM,  SUM,  CONCEN,  PI,  B,  STEP,  H 
INTEGER  XREF,  YREF,  PUFF,  TIM,  GRIDX,  GRIDY,  INCR 
C 

C  DIMENSION  STATEMENTS 

REAL  SIGY(60),  SIGX(60),  SIGZ(60),  A(60) ,  C(60),  D(60) 

REAL  DIST{60),  SOURCE{60) 

REAL  CONC(-101:101,  -101:101) 

REAL  WSPEED (60, 4) 

C 

C  CONSTANTS 

PI  =  3.141592654 
C 

C  INTRODUCTION  OF  THE  OVERALL  FRAME  OF  REFERENCE 

C  XREF=>  downwind  distance  (METERS)  from  common  frame  of  reference 
C  X  =>  XREF  *  100 

C  YREF=>  crosswind  distance  (METERS)  from  common  frame  of  reference 
C  Y  =>  YREF  *  100 

C  BETA=>  angle  between  reference  coordinate  and  the  reference  X-axis 
C  ALPHA=>  wind  direction  ("to-sees") 

C  The  Z  value  is  zero  for  ground- level  receptor  (METERS) 

DO  100  XREF  =  -GRIDX, GRIDX 
X  =  XREF  *  STEP 
C 

DO  200  YREF  =  GRIDY, -GRIDY, -1 
Y  =  YREF  *  STEP 
C 

Z  =  0. 

C 

SQ  =  (X**2)  +  (Y**2) 

HYP  =  SQRT  (SQ) 

C 

C  BETA  =>  ANGLE  BETWEEN  THE  REFERENCE  X-AXIS  AND  THE  RECEPTOR  POSITION 
C  Want  BETA  to  be  positive 

BETA  =  ATAN2D(Y,X) 

IF  (BETA.LT. 0 . 0)  BETA  =  BETA  +  360. 

C 

C  DETERMINATION  OF  WHETHER  IT  IS  NECESSARY  TO  CALCULATE  CONCENTRATION 
C  If  the  difference  between  the  wind  direction  and  the  angle  to  the 
C  receptor  falls  in  the  third  or  fourth  quadrant  of  the  shifted 

C  axis  the  concentration  will  be  set  to  zero 

GAMMA  =  ALPHA  -  BETA 
AGAMMA  =  ABS (GAMMA) 

IF  (GAMMA.LT. 0 . 0)  GAMMA  =  GAMMA  +  360. 

IF  ( (AGAMMA. GE. 90.) .AND. (AGAMMA. LE. 270.) )  GOTO  1000 

C 
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C  THE  SIGN  OF  THE  RECEPTOR  IS  NOT  IMPORTANT  HERE 
C  Depending  upon  the  quadrant  of  interest,  the  coordinates  of  the 
C  reference  frame  are  used  to  determine  the  coordinates  on  the 
C  shifted  frame  (based  upon  the  wind  direction) 

C  The  concentration  at  the  negative  YSHIFT  coordinates  are  the  same 
C  as  the  positive  YSHIFT  coordinates 
XSHIFT  =  HYP  *  COSD (GAMMA) 

YSHIFT  =  HYP  *  SIND (GAMMA) 

C 

C 

C  CALCULATION  OF  THE  CONCENTRATION  AT  RECEPTOR  POINT 
C 

C  CALCULATION  OF  THE  SIGMA  VALUES 

SIGY(PUFF)  =  A(PUFF) * (DIST(PUFF) **B) 

SIGX(PUFF)  =  SIGY(PUFF) 

SIGZ(PUFF)  =  C (PUFF) * (DIST(PUFF) **D (PUFF) ) 

C 

C 

C  NOW  FOR  THE  PUFF  EQUATION 

STRM  =  SOURCE (PUFF) /(( (2. *PI) ** (1.5) ) *SIGX (PUFF) * 

+  SIGY (PUFF) *SIGZ (PUFF) ) 

XTRM  =  -( (XSHIFT- (WSPEED(PUFF, 1) *TIM) ) **2) / 

+  (2  .* (SIGX(PUFF) **2) ) 

YTRM  =  - (YSHIFT**2)/(2.*(SIGY(PUFF)**2) ) 

ZTRM  =  - ( (H) **2) /(2.* (SIGZ(PUFF) **2) ) 

SUM  =  XTRM  +  YTRM  +  ZTRM 
IF(SUM.LT.-50.)  SUM  =  -50.0 
CONCEN  =  2* STRM* EXP (SUM) 

IF(CONCEN.LE. (O.lE-10) )  GOTO  1000 
GOTO  1100 
C 

C  DESIGNATION  OF  CONC  =  0  FOR  THE  APPROPRIATE  CONDITIONS 
1000  CONCEN  =0.0 

C 

C  CONCENTRATION  IS  CALCULATED  BY  ADDING  THE  AVERAGE  CONCENTRATION 
C  DURING  THE  MINUTE  TO  PREVIOUS  CALCULATED  CONCENTRATION 
1100  CONC(XREF,YREF)  =  CONC (XREF, YREF)  +  CONCEN/INCR 

C 

C  GO  BACK  TO  GET  Y  GRID  VALUES 
200  CONTINUE 

C  GO  BACK  TO  GET  X  GRID  VALUES 
100  CONTINUE 
C 
C 

RETURN 

END 
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SUBROUTINE  RPTDATA  (STN,  YR,  DATE,  TIME,  SP,  A,  DTDZ,  ZSP,  HGT, 

+  TEMP,  RH,  WSPEED,  WDIR,  STHETA,  P) 

C 

C  THIS  SUBROUTINE  READS  IN  PUFF  DATA  THAT  IS  THE  SAME  AS  THE  PREVIOUS 
C 

C  DECLARATION  OF  VARIABLE  TYPES 
INTEGER  P 
C 

C  DIMENSION  STATEMENTS 

REAL  TEMP(60,4),  RH(60,4),  WSPEED (60, 4) ,  WDIR(60,4) 

REAL  STHETA (6 0,4) 

REAL  SP(60),  A(60),  DTDZ(60),  ZSP{60) 

INTEGER  STN(60),  YR{60),  DATE(60),  TIME(60) 

INTEGER  HGT (60, 4) 

C 

STN(P)  =  STN(P-l) 

YR(P)  =  YR{P-1) 

DATE(P)  =  DATE (P-1) 

TIME (P)  =  TIME (P-1) 

SP(P)  =  SP(P-l) 

A(P)  =  A(P-l) 

DTDZ (P)  =  DTDZ (P-1) 

ZSP(P)  =  ZSP(P-l) 

DO  100  J  =  1,4 

HGT(P,J)  =  HGT(P-1,J) 

TEMP(P,J)  =  TEMP(P-1,J) 

RH(P,J)  =  RH(P-1,J) 

WSPEED(P,J)  =  WSPEED(P-1, J) 

WDIR(P,J)  =  WDIR(P-1,J) 

STHETA(P,J)  =  STHETA(P-1, J) 

100  CONTINUE 

RETURN 
END 


A-22 


SUBROUTINE  STABPAR  (STHETA,  TEMP,  ZSP,  P,  A,  B) 

C 

C  THIS  SUBROUTINE  CALCULATES  VERTICAL  AND  HORIZONTAL  STABILITY 
C  PARAMETERS  AND  HORIZONTAL  COEFFICIENT  AND  EXPONENT 
C 

C  DECLARATION  OF  VARIABLE  TYPES 
INTEGER  P 

REAL  TEMPI,  TEMPIO,  DT,  DZ,  DTDZ,  B 
C 

C  DIMENSION  STATEMENTS 

REAL  SP{60),  ZSP{60),  A(60) 

REAL  STHETA(60,4) ,  TEMP(60,4) 

C 

IF  {STHETA{P,1)  .GT.  (7.5) )  THEN 

SP(P)  =  (-0.2)  *  STHETA(P,1)  +  5.5 
ELSEIF  (STHETA(P, 1)  .EQ.  (7.5) )  THEN 
SP(P)  =4.0 

ELSEIF  (STHETA(P,1) .GT. (3.8) )  THEN 
SP(P)  =  (-0.27  *  STHETA(P,1))  +  6.03 
ELSEIF  (STHETA(P,1)  .EQ.  (3.8) )  THEN 
SP(P)  =5.0 

ELSEIF  (STHETA(P,1)  .GT.  (2.95) )  THEN 
SP(P)  =  (-0.59  *  STHETA(P,1))  +  7.23 
ELSEIF  (STHETA(P, 1)  .LE.  (2.95) )  THEN 
SP(P)  =  5.5 
END  IF 

C  CALCULATION  OF  THE  VERTICAL  STABILITY  PARAMETER 
C  Conversion  of  temperature  from  Fahrenheit  to  Celcius 
C  NRC  Method  of  classification  is  based  upon  DT/DZ  (C/lOOm) ; 

C  therefore,  need  to  convert  by  multiplying  temp  by  (9/100) 

C 

TEMPIO  =  (5./9.)*  TEMP(P,1)  -  32. 

TEMPI  =  (5./9.)*  TEMP(P,4)  -  32. 

DT  =  (TEMPIO  -  TEMPI) 

DZ  =  9.0 

DTDZ  =  (9,/100.)* (DT/DZ) 

C 

IF  (DTDZ. LT. -1.9)  THEN 
ZSP(P)  =  0.5 

ELSEIF  (DTDZ. LE. -1.5)  THEN 
ZSP(P)  =  5*DTDZ  +  10.5 
ELSEIF  (DTDZ. LE. -0.5)  THEN 
ZSP(P)  =  DTDZ  +4.5 
ELSEIF  (DTDZ. LE. 1.5)  THEN 
ZSP(P)  =  0.5*DTDZ  +  4.25 
ELSEIF  (DTDZ. GE. 1.5)  THEN 
ZSP(P)  =  5.5 
END  IF 

A(P)  =  0.479  -  (0.1231*SP(P) )  +  (0 . 00904* (SP (P) **2 .) ) 

B  =  0.9 

RETURN 

END 
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SUBROUTINE  TOPTEN(P,  GRIDX,  GRIDY,  CONG,  TOPIO,  COORDX,  COORDY, 
+  WHEN) 

C  THIS  SUBROUTINE  CALCULATES  THE  "TOP  10"  CONCENTRATIONS 
C 

C  DECLARATION  OF  VARIABLE  TYPES 

INTEGER  GRIDX,  GRIDY,  P,  TEMPX,  TEMPY 
REAL  TEMP 
C 

C  DIMENSION  STATEMENTS 

REAL  CONC(- 101: 101,  -101:101),  TOPIO(IO) 

INTEGER  COORDX  (10),  COORDY (10),  WHEN (10) 

C 

DO  200  YREF  =  GRIDY,  -GRIDY,  -1 
DO  300  XREF  =  -GRIDX,  GRIDX 
C 

C  IF  VALUE  IS  LESS  THAN  TENTH  PLACE,  GO  TO  THE  NEXT  VALUE 
IF  (CONC (XREF, YREF) .LE. TOPIO (10) ) GOTO  300 

C  THIS  LOOP  DETERMINES  THE  POSITION  OF  THE  NEW  VALUE 
DO  400  I  =  1,10 

IF  (CONC (XREF, YREF) .GT.TOPIO (I) )  GOTO  425 
400  CONTINUE 

GOTO  300 


C 


NOW  PULL  DOWN  APPROPRIATE  NUMBER  OF  VALDES  TO  MAKE  ROOM  FOR  NEW 
425  DO  500  J=9,I,-1 

TOPIO (J+1)  =  TOPIO (J) 

COORDX ( J+1 )  =  COORDX ( J) 

COORDY (J+1)  =  COORDY (J) 


WHEN (J+1) 

=  WHEN(J) 

CONTINUE 

TOPIO  (I)  = 

CONC (XREF, YREF) 

COORDX (I)  = 

XREF 

COORDY ( I )  = 

YREF 

WHEN (I) 

P 

C 


300  CONTINUE 
200  CONTINUE 


C 


RETURN 

END 
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Appendix  B  *METALU  FORTRAN  77  Code 


program  metal 1 

integer  yr, time, tper, hgt (4) 

C  DIMENSION  HEIGHT,  ETC--FOUR;  WX  VARIABLES --FOUR  BY  SIXTY 

dimension  sigth2 (4) ,  wddif2 (4) ,  avexcom(4) , 
aveycom(4) ,  rhtot(4) 

dimension  avetemp(4),  averh{4),  avewd(4) ,  avews(4), 
avesig (4) 

dimension  ycomp(4),  xcomp(4),  xtot(4),  ytot(4), 
wstot(4),  temptot(4) 

dimension  temp (4, 60),  relhum(4, 60) ,  ws(4,60), 
wd(4,60),  stheta{4,60) 
dimension  oldwd(4) 

open  (unit=7,  f ile= ' s5061001 .dat ' ,  status= ' old' ) 
open  (unit =8,  file= ' average ' ,  status= 'unknown' ) 

open  (unit=9,  f ile= ' s560min' ,  status= 'unknown' ) 

C  SELECT  THE  TIME  PERIOD  OVER  WHICH  THE  AVERAGE  WILL  BE 
C  TAKEN:  TPER  SET  THE  COUNTERS  WSTOT,  XTOT,  AND  YTOT  TO  ZERO 

tper  =  60 

do  55  i=l,4 
wstot(i)  =  0.0 
xtot(i)  =  0.0 
rhtot(i)  =  0.0 
ytot{i)  =  0.0 
55  continue 

c  COLLECT  ALL  THE  WEATHER  DATA  FOR  THE  PERIOD 
1  do  100  j=l,tper 

read(7, *,end=999) id,  yr,  jd,  time,  (hgt(i),  temp(i,j), 

relhum(i,j),  ws(i,j),  wd(i,j), 
stheta{i,j),  i=l,4) 
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C  PERFORM  SOME  OF  THE  PRELIMINARY  CALCULATIONS 
c  BEGIN  BY  OBTAINING  THE  X  AND  Y  COMPONENTS  OF  THE  WIND 
c  DIRECTION 

c  MULTIPLY  BY  0.0174533  TO  CONVERT  DEGREES  TO  RADIANS  FOR 
C  TRIG  FUNCTIONS 

do  45  i=l,4 

xcomp(i)  =  cos ( (wd(i, j ) -270 . 0) *0 . 0174533) 

ycoinp(i)  =  -1.0  *  sin  (  (wd  (i,  j  ) -270 . 0)  *0 . 0174533) 

C  TOTAL  THE  X  AND  Y  COMPONENTS  FOR  WIND  DIRECTION  AND 
C  WSTOT  FOR  WIND  SPEED,  TEMPTOT  FOR  TEMPERATURE,  AND  RHTOT 
C  FOR  REL  HUMIDITY 

xtot(i)  =  xtot(i)  +  xcomp{i) 
ytot(i)  =  ytot(i)  +  ycomp(i) 
wstot(i)  =  wstot(i)  +  ws(i,j) 
temptot{i)  =  temptot(i)  +  temp(i,j) 
rhtot(i)  =  rhtot(i)  +  relhum(i,j) 

C  RETURN  TO  READ  NEXT  LINE  OF  DATA 

45  continue 

C  WRITE  THE  PRELIMINARY  OUTPUT  FOR  LATER  COMPARISON 

write (8 , 10) jd,  time, (hgt (i) ,  temp(i,j),  relhum(i,j), 
ws(i,j),  wd(i,j),  stheta(i,j),  i=l,4) 

10  format(I3,  14,  4(12,  3F5.1,  F6.1,  F5.1)) 

100  continue 

C  HAVING  COMPLETED  THE  LOOP,  FOR  EACH  LEVEL 
c  CALCULATE  FINAL  VALUES  FOR  X  AND  Y  COMPONENTS  OF  WIND 
C  DIRECTION,  THE  AVERAGE  WIND  SPEED,  AVERAGE  REL  HUM,  AND 
c  THE  AVERAGE  TEMPERATURE 

c  ALSO  IN  THAT  LOOP,  RESET  RHTOT,  TEMPTOT,  WSTOT,  XTOT,  AND 
C  YTOT  TO  ZERO 

do  65  i=l,4 
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avexcom(i)  =  xtot{i)  /  tper 
aveycom(i)  =  ytot(i)  /  tper 
avews(i)  =  wstot(i)  /  tper 
avetemp(i)  =  temptot(i)  /  tper 
averh(i)  =  rhtot(i)  /  tper 
if (averh(i) .ge . 100 . 0)  averh(i)  =  99.9 

xtot(i)  =  0.0 
ytot(i)  =  0.0 
wstot(i)  =  0.0 
temptot(i)  =  0.0 
rhtot(i)  =  0.0 

c  DETERMINE  THE  AVERAGE  WIND  DIRECTION  FROM  XFIN  AND  YFIN 
C  REVERSE  ARGUEMENTS  TO  CONFORM  WITH  CARTESIAN  COORDINATES 
c  THEN  ADD  180  DEGREES  TO  ACCOUNT  FOR  FACT  THAT  WINDS  ARE 
c  "FROMSEES" 

avewd{i)  =  atan2d  (avexcom(i)  ,  aveycoin(i)  )  +  180.0 

C  THEN  GET  RID  OF  NEGATIVE  WIND  DIRECTIONS 

if (avewd (i) . It . 0 . 0)  avewd(i)  =  avewd(i)  +  360.0 

C  IF  WIND  SPEED  IS  ZERO,  AVERAGE  WIND  DIRECTION  IS  CARRIED 
C  FROM  PREVIOUS 

if(avews(i)  .eq.  0.0)  avewd (i)  =  oldwd(i) 

C  CALCULATE  THE  AVERAGE  SIGMA  THETA  FOR  PERIOD 
C  FIRST  INITIALIZE  COUNTERS  SIG2TH  AND  WDDIF2 

sigth2 (i)  =0.0 
wddif2{i)  =  0.0 

C  SECOND:  SUM  THE  SQUARES  OF  THE  INDIVIDUAL  SIGMA  THETA 
C  VALUES 

do  25  j  =  l,tper 

sigth2(i)  =  sigth2 (i)  +  stheta(i,j)  *  stheta(i,j) 
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C  NEXT  DETERMINE  THE  DIFFERENCE  BETWEEN  INDIVIDUAL  WIND 
C  DIRECTIONS  AND  MEAN  WIND  DIRECTION  FOR  PERIOD.  REDUCE  ALL 
c  VALUES  GREATER  THAN  180  DEGREES,  THEN  USE  SUPPLEMENTAL 
C  ANGLE.  THEN  SUM  SQUARES. 

If{ws(i,j)  .eq.  0.0)  wd(i,j)  =  avewd{i) 
wd ( i , j )  =  avewd ( i )  -  wd ( i , j ) 

if(wd(i,j)  .It.  -180.0)  wd(i,j)  =  360.0  +  wd(i,j) 
if(wd(i,j)  .gt.  180.0)  wd(i,j)  =  360.0  -  wd(i,j) 
wddif2(i)  =  wddif2(i)  +  wd(i,j)  *  wd(i,j) 

25  continue 

C  FINALLY  TAKE  SUM  OF  SQUARES  OF  SIGMAS  AND  WIND 
C  DIFFERENCES,  DIVIDE  BY  NUMBER  OF  OBSERVATIONS  TO  GET  THE 
C  NEW  SIGMA  THETA. 

avesig(i)  =  sqrt ( (sigth2 (i)  +  wddif2(i))  /  tper) 
oldwd(i)  =  avewd (i) 

65  continue 

C  WRITE  OUTPUT  FOR  THIS  PERIOD 

write (8 , 70) jd, time,  (hgt (i) , avetemp (i) , averh (i)  , 

&:avews  (i)  ,  avewd  (i)  ,avesig{i)  ,  1=1,4) 

70  format(i4,i4,4(I2,  3F5.1,  F6.1,  F5.1)) 


write (9 , 75) id, yr , jd, time, (hgt (i) , avetemp (i) , averh (i) , 
&avews (i) , avewd(i) , avesig(i)  ,  1=1, 4) 

75  format (i2, i5, i4, i5, lx,4 (12,  3F5.1,  F6.1,  F5.1)) 

C  RETURN  TO  BEGINNING  AND  PROCEED  WITH  NEXT  AVERAGING  PERIOD 

go  to  1 

999  Stop 
end 
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Appendix  C  Derivation  of  Calculation  Used  to  Determine  Sigma  Theta 


a0  (horizontal  standard  deviation  of  the  wind)  for  each  minute  of  data  gathered  from 
meteorological  towers 

Oct 

ae  for  averaging  time  of  interest,  i.e.,  10  or  60  minutes 

X 

individual  one  second  wind  direction  readings  taken  by  computer  in  meteorological 
towers  (one  second  readings  are  not  available  as  output  data) 

XT 

individual  one-minute  wind  direction  readings  taken  by  computer  in  meteorological 
towers 

mean  one-minute  wind  direction  calculated  by  meteorological  towers 

H 

grand  mean  of  10  or  60  minute  wind  direction 

N 

number  of  one  second  readings  (60)  used  to  calculate  the  mean  wind  direction 

Nt 

number  of  one  second  readings  needed  to  calculate  the  10  or  60  minute  averages 
(600  readings  for  10  minutes  -and-  3600  readings  for  60  minutes) 

Nr 

ratio  N/Nx 

Calculation  of  Standard  Deviations 

For  one  minute  readings 

a 


^1 


0-1 


^  2 
^1 


2  _  E  F\y 


N 


z 


N 


N 


Ml 


Z =  {-I  ^  m^i)n 


For  1 0  and  60  min  ute  readings 
Calculation  of  the  "  grand" 
standard  deviation 


2 

Nj- 

_2  Z(^T-2xTMg+Ml) 

^  cr 


= 


2 

T 


c-i 


From  the  above  example  for  ^ 

Substituting  and  remembering Xj  =  ^x 

iV7’(o-|  +  /i|)  =  2  (o-f  +//f)v 


'y 

Simplifying  by  solving  fora  „ 


Taking  the  square  root  of  both  sides  and  expanding 


Since 


This  equation  for  Cg  was  used  to  calculate  the  o  for  the  appropriate  averaging  time,  i.e.,  1, 
10,  or  60  minutes,  in  the  FORTRAN  77  code  for  ‘METALL’  (See  Appendix  B). 

*See  Proof  on  the  next  page 
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♦Based  upon  “  <“|) 

Expanding 


A 


-  2  p  gY,  p  +  Y,  mI  =  ^ 

Sim plifying  and  combining  terms 

=  -2E  g 

Simplifying 

Z  =  Z  f^g 

M  =  ffM  g 

therefore , 
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